; MIDI_CHPDAT - Reduce chopped data from MIDI. ; chpres = midi_chpdat(rawdat, err_stat) ; ; MIDI_WAVAXI - Plot a wavelength axis on top of the previous plot. ; midi_wavaxi, disper, xtickl=0.02, /pltlab ; ; MIDI_OUTDIS - Plot one detector window for dispersed MIDI data. ; midi_outdis, detwin, positi, offset, ctable=3, disper='prism' ; ; MIDI_OUTFIE - Plot one detector window for undispersed MIDI data. ; midi_outfie, detwin, positi, offset, ctable=3, /llabel, /rlabel ; ; MIDI_PLTWIN - Make a plot of the two detector windows of MIDI data. ; midi_pltwin, rawdat, psfile=psfile, /sqroot, /invers, factor=3 ; ; MIDI_PHOTOM - Chop a MIDI photometry and extract the spectrum. ; photom = midi_photom(phtdat, /output, /nosky, phtbnd=[13,21]) ; ; MIDI_PHTFIL - Reduce a MIDI photometry data file. ; photom = midi_phtfil(phtfil, /nosky, phtbnd=[13,21], skydis=3) ; ; MIDI_PLTFRI - Show a video of the delay line position and the fringe. ; midi_pltfri, ifdata, frames, delay, tsmoot=20, /davera, factor=3 ; ; MIDI_GAUDIS - Fit a Gaussian to every channel in a dispersed MIDI window. ; fitpar = midi_gaudis(detwin, rrange=[9,26], /output, /quiet) ; ; MIDI_MSKFIT - Create a MIDI mask from A and B photometry files. ; mskdat = midi_mskfit(infile, msknam, factor=1.25, /output) ; ; MIDI_FLAGOK - Determine which frames where flagged as good in a fringe track. ; flagok = midi_flagok(flgdat, gdelay) ; ; MIDI_UNWRAP - Unwrap phases assuming a smooth phase change ; newphi = midi_unwrap(oldphi, smooth, maxval=40) ; ; MIDI_GETPHI - Get the water vapour phases from a FITS file. ; phases = midi_getphi('filnam', /unwrap) ; ; MIDI_PLTPHI - Plot the water vapour phases for a MIDI fringe track. ; midi_pltphi, filnam, xrange=[0,9], yrange=[0,1], psfile='psfile' ; ; MIDI_PLTOPD - Plot the Optical Path Difference for a MIDI fringe track. ; midi_pltopd, filnam, xrange=[0,9], yrange=[0,1], psfile='psfile' ; ; MIDI_PLTOVL - Plot the beam overlap for the photometry, fringes and masks. ; midi_pltovl, datset, [11,13], /normal, psfile='/tmp/overlap' ; MIDI_PLTRES - ; ; MIDI_PLTRAD - Plot the correlated flux or visibility vs. the baseline length. ; midi_pltrad, caldat, 12., 0.2, /visamp, xrange=[0,90], /ylog ; ; MIDI_PLTUVP - Plot points in a UV plane with values from a MIDI result struct. ; midi_pltuvp, caldat, wavele, tagnam='corflx', zrange=[0,10] ;+ ; NAME: ; MIDI_CHPDAT ; ; PURPOSE: ; Reduce chopped data from MIDI. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; chpres=midi_chpdat(rawdat [, err_stat]) ; ; REQUIRED INPUTS: ; rawdat structure containing the chopped data, e.g. the photometry ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; none ; ; OUTPUTS: ; chpres structure containing the difference of object - sky ; err_stat ; ; DESCRIPTION AND EXAMPLE: ; The programme has to be provided with a structure RAWDAT containing ; chopped data from MIDI. Both dispersed data (e.g. photometry) or ; undispersed data (e.g. acquisition) is accepted. For example chopped ; data can be loaded with oirGetData: ; rawdat = oirGetData('MIDI.2000-01-01T12:34:56.000.fits'). ; The programme then determines which frames are sky frames and which ; frames are target frames from the tag TARTYP in the structure. ; Then the mean of the sky frames is subtracted from the mean of the ; target frames for both detector windows. The result is returned by ; the programme as a structure of the same format as the input data, ; although the data fields are converted to float if necessary. The ; simplest way to run the programme is by calling; ; chpres=midi_chpdat(rawdat) ; ; CALLED BY: ; none ; ; CALLING: ; mod_struct ; ; MODIFICATION HISTORY: ; 2008-03-15 Written by Konrad R. Tristram ; 2010-11-03 Konrad R. Tristram: Rewrite to cope with multiple windows. ; FUNCTION MIDI_CHPDAT, RAWDAT, ERR_STAT ; CHECK IF SKY AND TARGET FLAGS HAVE BEEN SET CORRECTLY ;------------------------------------------------------------------------------- catch, err_stat if err_stat ne 0 then begin message, 'Something is wrong with the data. Returning.', /continue return, -1 endif ; FIND HOW MANY AND WHERE THE DATA WINDOWS ARE ;------------------------------------------------------------------------------- tagnam = tag_names(rawdat) winidx = where(strmatch(tag_names(rawdat), 'DATA*')) if winidx[0] lt 0 then begin message, 'No detector data fields detected in the data. Exiting!', /cont return, -1 endif winnum = n_elements(winidx) ; USE THE FIRST FRAME AS THE RESULT STRUCT ;------------------------------------------------------------------------------- result = rawdat[0] ; LOOP OVER THE WINDOWS ;------------------------------------------------------------------------------- for i=0,winnum-1 do begin ; FIND WHERE THE TARGET TYPE FOR THIS WINDOW IS LOCATED ;----------------------------------------------------------------------- ; The follwing line is a workaround for SCI_PHOT data if winnum eq 4 then tmpvar = (i+1)>2<3 else tmpvar = i+1 taridx = where(tag_names(rawdat) eq 'TARTYP'+strtrim(tmpvar, 2)) ; FIND THE FRAMES WHICH ARE ON TARGET AND THOSE WHICH ARE ON SKY ;----------------------------------------------------------------------- tar = where(rawdat.(taridx) eq 'T ') sky = where(rawdat.(taridx) eq 'S ') ; CALCULATE THE TOTALS AND SUBTRACT THE SKY FROM THE OBJECT ;----------------------------------------------------------------------- tmpvar = total(rawdat[tar].(winidx[i]), 3) / float(n_elements(tar)) - $ total(rawdat[sky].(winidx[i]), 3) / float(n_elements(sky)) ; SAVE THE DATA IN THE RESULT STRUCT; WHICH HAS TO BE CHANGED TO FLOAT ;----------------------------------------------------------------------- result = mod_struct(result, tagnam[winidx[i]], tmpvar) endfor ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, result END ;+ ; NAME: ; MIDI_WAVAXI ; ; PURPOSE: ; Plot a wavelength axis on top of the previous plot. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; midi_wavaxi, disper [, xtickl=xtickl] [, /pltlab] ; ; REQUIRED INPUTS: ; disper a string specifying the dispersive element ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; xtickl length of the tick marks ; pltlab also plot tick labels ; ; OUTPUTS: ; a plot of an alternative x-axis with wavelength marks ; ; DESCRIPTION AND EXAMPLE: ; This programme plots an alternative x-axis to a previously generated ; plot. The previous plot has to have an x-axis with the pixel positions ; on the MIDI detector when read out in dispersed mode. The new axis has ; major tick marks at wavelengths from 6 to 14 micron in 1 micron steps ; and minor tickmarks in 0.25 micron steps. The position of the tick ; marks is determined by a second order polynamial which depends on ; DISPER. The coefficients for the polinomial are defined for ; DISPER='prism' and DISPER='grism'. The values of the coefficients are ; taken from the MIDI manual (Issue 82, Date 29/02/2008). If the keyword ; /PLTLAB is set, also the values of the major tickmarks are plotted. ; The length of the tickmarks as a fraction of the flot hight can be ; specified with XTICKL. The simplest way to run the programme is by ; calling: ; midi_wavaxi, 'prism' ; ; CALLED BY: ; midi_outdis ; midi_photom ; midi_mskfit ; midi_pltres ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2008-03-15 Written by Konrad R. Tristram ; PRO MIDI_WAVAXI, DISPER, XTICKL=XTICKL, PLTLAB=PLTLAB if not keyword_set(xtickl) then xtickl=0.02 ; DEFINE THE COEFFICIENTS OF THE WAVELENGTH CONVERSION (FROM THE MIDI MANUAL) ;------------------------------------------------------------------------------- if strcmp('prism', disper, 5, /fold) then begin a_wpol = -1.539e-4 b_wpol = 9.491e-3 c_wpol = 15.4519 d_wpol = 120.0000 s_wpol = -1.0000 endif else if strcmp('grism', disper, 5, /fold) then begin a_wpol = -1.211e-6 b_wpol = 23.122e-3 c_wpol = 7.5442 d_wpol = 5.0000 s_wpol = +1.0000 endif else begin a_wpol = 0.0000 b_wpol = 1.0000 c_wpol = 0.0000 d_wpol = 0.0000 s_wpol = +1.0000 endelse ;waveLengthArray = oirGetWaveLength(detectorFile) ; PLOT MINOR TICKS FOR THE AXIS IN 0.25 MICRON INTERVALS ;------------------------------------------------------------------------------- waveletk = indgen(37)/4.+5. waveletk = waveletk[where(waveletk mod 1 ne 0)] waveleps = (-b_wpol + s_wpol * sqrt(b_wpol^2-4.*a_wpol*(c_wpol-waveletk))) / $ 2. / a_wpol - d_wpol axis, xaxis=1, xstyle=1, xtickv=waveleps, xticklen=0.5*xtickl, $ xtickname=replicate(' ', n_elements(waveletk)), $ xticks=n_elements(waveletk)-1 ; PLOT MAIN TICKS AND TICK LABELS (IF KEYWORD SET) IN 1 MICRON INTERVALS ;------------------------------------------------------------------------------- waveletk = indgen(10)/1.+5. waveleps = (-b_wpol + s_wpol * sqrt(b_wpol^2-4.*a_wpol*(c_wpol-waveletk))) / $ 2. / a_wpol - d_wpol if keyword_set(pltlab) then begin axis, xaxis=1, xstyle=1, xtickv=waveleps, xticklen=xtickl, $ xtickname=strcompress(string(waveletk, format='(I2)'), /remove), $ xticks=n_elements(waveletk)-1 endif else begin axis, xaxis=1, xstyle=1, xtickv=waveleps, xticklen=xtickl, $ xtickname=replicate(' ', n_elements(waveletk)), $ xticks=n_elements(waveletk)-1 endelse END ;+ ; NAME: ; MIDI_OUTDIS ; ; PURPOSE: ; Make a plot of one detector window for dispersed MIDI data. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; midi_outdis, detwin, positi [, offset] [, ctable=ctable] ; [, disper=disper] [, /tlabel] [, /blabel] [, /llabel] ; [, /rlabel] ; ; REQUIRED INPUTS: ; detwin a single dispersed MIDI detector window ; positi the position of the plot in device coordinates ; ; OPTIONAL INPUTS: ; offset the pixel position for zero spatial offset ; ; KEYWORDS: ; ctable the colour table to be used for plotting the array ; disper the dispersion type, e.g. 'prism' or 'grism' ; tlabel plot tick labels for x-axis on top of the plot ; blabel plot tick labels for x-axis on bottom of the plot ; llabel plot tick labels for left y-axis ; rlabel plot tick labels for right y-axis ; ; OUTPUTS: ; a plot of the provided array surrounded by multiple axes ; ; DESCRIPTION AND EXAMPLE: ; This programme creates a plot for the provided single MIDI detector ; window for dispersed data (e.g. photometry). The detector array must be ; previously scaled correclty, that is 0 < DETWIN < 255. The position of ; the plot in device coordinates is governed by POSITION. The x- and ; y-axes are scaled with the pixel position on the detector; the ; alternative x-axis is scaled with the wavelength and depends on the ; keyword DISPER; the alternative y-axis is scaled in arcseconds with the ; pixel position of the zeropoint defined by OFFSET. For prism data, the ; simplest way to run the programme is by calling for example: ; midi_outdis, detwin, [50, 50, 3*171+50, 3*41+50], disper='prism' ; ; CALLED BY: ; midi_pltwin ; midi_photom ; midi_pltres ; ; CALLING: ; midi_wavaxi ; ; MODIFICATION HISTORY: ; 2008-03-15 Written by Konrad R. Tristram ; 2009-06-30 Konrad R. Tristram: Added keywords LLABEL and RLABEL. ; PRO MIDI_OUTDIS, DETWIN, POSITI, OFFSET, CTABLE=CTABLE, DISPER=DISPER, $ TLABEL=TLABEL, BLABEL=BLABEL, LLABEL=LLABEL, RLABEL=RLABEL ; SET DEFAULT VALUES FOR VARIABLES IF NOT PROVIDED BY THE USER ;------------------------------------------------------------------------------- if n_elements(offset) ne 1 then offset=14.5 if n_elements(ctable) eq 0 then ctable = 3 $ else ctable = (byte(ctable))[0] < 45 if n_elements(disper) ne 1 then disper='unknown' ; SET THE PIXEL SCALE AND GET THE ARRAY SIZE ;------------------------------------------------------------------------------- pscale = 0.086 ; from Tristram 2007, PHD-Thesis arrsiz = size(detwin) ; DETERMINE IF THE ARRAY HAS TO BE MAGNIFIED AND CALCULATE THE NEW SIZES ;------------------------------------------------------------------------------- if !d.name eq 'X' or !d.name eq 'WIN' then begin x_scal = positi[2]-positi[0] y_scal = positi[3]-positi[1] endif else begin x_scal = n_elements(detwin[*,0]) y_scal = n_elements(detwin[0,*]) endelse ; PLOT OUT THE ARRAY USING THE SPECIFIED COLOUR TABLE ;------------------------------------------------------------------------------- loadct, ctable, /silent tv, congrid(detwin, x_scal, y_scal), positi[0], positi[1], $ xsize=positi[2]-positi[0], ysize=positi[3]-positi[1] loadct, 39, /silent ; ADD THE PIXEL AXES AND CREATE THE BOTTOM LABEL (IF KEYWORD IS SET) ;------------------------------------------------------------------------------- if keyword_set(llabel) then ytickn = '' else ytickn = replicate(' ', 50) if keyword_set(blabel) then begin plot, [0,0], [0,0], xrange=[-0.5, arrsiz[1]-0.5], $ yrange=[-0.5, arrsiz[2]-0.5], xstyle=9, ystyle=9, $ xtickv=waveleps, xtitle='pixels on detector', $ ytickname=ytickn, xticklen=0.08, yticklen=0.02, yminor=5, $ position=positi, /device, /nodata, /noerase endif else begin plot, [0,0], [0,0], xrange=[-0.5, arrsiz[1]-0.5], $ yrange=[-0.5, arrsiz[2]-0.5], xstyle=9, ystyle=9, $ xtickv=waveleps, xtickname=replicate(' ', 50), $ ytickname=ytickn, xticklen=0.08, yticklen=0.02, yminor=5, $ position=positi, /device, /nodata, /noerase endelse ; ADD THE WAVELENGTH AXES ON TOP OF THE PLOT ;------------------------------------------------------------------------------- if keyword_set(tlabel) then midi_wavaxi, disper, xtickl=0.08, /pltlab $ else midi_wavaxi, disper, xtickl=0.08 ; PLOT THE RIGHT Y-AXIS WHICH IS SCALED WITH ANGULAR SCALE ;------------------------------------------------------------------------------- if keyword_set(rlabel) then ytickn = '' else ytickn = replicate(' ', 50) yrange = ([-0.5, arrsiz[2]-0.5] - offset) * pscale axis, yaxis=1, yrange=yrange, ystyle=1, ytickname=ytickn, yticklen=0.02, $ ytickin=1, yminor=5 END ;+ ; NAME: ; MIDI_OUTFIE ; ; PURPOSE: ; Make a plot of one detector window for undispersed MIDI data. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; midi_outfie, detwin, positi [, offset] [, ctable=ctable] ; [, /llabel] [, /rlabel] [, /tlabel] ; ; REQUIRED INPUTS: ; detwin a single undispersed MIDI detector window ; positi the position of the plot in device coordinates ; ; OPTIONAL INPUTS: ; offset the pixel positions for zero angular offsets in x and y ; ; KEYWORDS: ; ctable the colour table to be used for plotting the array ; llabel plot tick labels for left y-axis ; rlabel plot tick labels for right y-axis ; tlabel plot axis label for top x-axis ; ; OUTPUTS: ; a plot of the provided array surrounded by pixel and angular offset axes ; ; DESCRIPTION AND EXAMPLE: ; This programme creates a plot for the provided single MIDI detector ; window for undispersed data (e.g. acquisition). The detector array must ; be previously scaled correclty, that is 0 < DETWIN < 255. The position ; of the plot in device coordinates is governed by POSITION. The x- and ; y-axes are scaled with the pixel position on the detector; the ; alternative x- and y-axes are scaled in arcseconds with the pixel ; position of the zeropoints defined by the two element array given in ; OFFSET. For acquisition data, the simplest way to run the programme is ; by calling for example: ; midi_outfie, detwin, [50, 50, 3*62+50, 3*69+50], disper='prism' ; ; CALLED BY: ; midi_pltwin ; ; MODIFICATION HISTORY: ; 2008-03-15 Written by Konrad R. Tristram ; 2010-12-06 Konrad R. Tristram: Added keyword TLABEL. ; PRO MIDI_OUTFIE, DETWIN, POSITI, OFFSET, CTABLE=CTABLE, $ LLABEL=LLABEL, RLABEL=RLABEL, TLABEL=TLABEL ; SET DEFAULT VALUES FOR VARIABLES IF NOT PROVIDED BY THE USER ;------------------------------------------------------------------------------- if n_elements(offset) ne 2 then offset=[29.0,31.0] if n_elements(ctable) eq 0 then ctable = 3 $ else ctable = (byte(ctable))[0] < 45 ; SET THE PIXEL SCALE AND GET THE ARRAY SIZE ;------------------------------------------------------------------------------- pscale = 0.086 ; from Tristram 2007, PHD-Thesis arrsiz = size(detwin) ; DETERMINE IF THE ARRAY HAS TO BE MAGNIFIED AND CALCULATE THE NEW SIZES ;------------------------------------------------------------------------------- if !d.name eq 'X' or !d.name eq 'WIN' then begin x_scal = positi[2]-positi[0] y_scal = positi[3]-positi[1] endif else begin x_scal = n_elements(detwin[*,0]) y_scal = n_elements(detwin[0,*]) endelse ; PLOT OUT THE ARRAY USING THE SPECIFIED COLOUR TABLE ;------------------------------------------------------------------------------- loadct, ctable, /silent tv, congrid(detwin, x_scal, y_scal), positi[0], positi[1], $ xsize=positi[2]-positi[0], ysize=positi[3]-positi[1] loadct, 39, /silent ; ADD THE PIXEL AXES AND CREATE THE LEFT LABELS (IF KEYWORD SET) ;------------------------------------------------------------------------------- if keyword_set(llabel) then begin plot, [0,0], [0,0], xrange=[-0.5, arrsiz[1]-0.5], $ yrange=[-0.5, arrsiz[2]-0.5], xstyle=9, ystyle=9, $ ytitle='pixels on detector', xtitle='pixels on detector', $ xticklen=0.03, yticklen=0.03, xminor=5, yminor=5, $ xtickv=waveleps, position=positi, /device, /nodata, /noerase endif else begin plot, [0,0], [0,0], xrange=[-0.5, arrsiz[1]-0.5], $ yrange=[-0.5, arrsiz[2]-0.5], xstyle=9, ystyle=9, $ ytickname=replicate(' ', 50), xtitle='pixels on detector', $ xticklen=0.03, yticklen=0.03, xminor=5, yminor=5, $ xtickv=waveleps, position=positi, /device, /nodata, /noerase endelse ; PLOT TOP X-AXIS WITH ANGULAR SCALE ;------------------------------------------------------------------------------- xrange = ([-0.5, arrsiz[1]-0.5] - offset[0]) * pscale axis, xaxis=1, xrange=xrange, xstyle=1, xticklen=0.03, $ xtickin=1, xminor=5 ; PLOT THE LABEL FOR THE TOP X-AXIS ;------------------------------------------------------------------------------- if keyword_set(tlabel) then xyouts, (positi[2]+positi[0])/2, $ 0.5 * positi[1] + positi[3], 'offset [arcsec]', align=0.5, /device ; PLOT RIGHT Y-AXIS AND CREATE THE RIGHT LABELS (IF KEYWORD SET) ;------------------------------------------------------------------------------- yrange = ([-0.5, arrsiz[2]-0.5] - offset[1]) * pscale if keyword_set(rlabel) then begin axis, yaxis=1, yrange=yrange, ystyle=1, yticklen=0.03, $ ytitle='offset [arcsec]', ytickin=1, yminor=5 endif else begin axis, yaxis=1, yrange=yrange, ystyle=1, yticklen=0.03, $ ytickname=replicate(' ', 50), ytickin=1, yminor=5 endelse END ;+ ; NAME: ; MIDI_PLTWIN ; ; PURPOSE: ; Make a plot of the detector windows of MIDI data. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; midi_pltwin, rawdat [, ctable] [, psfile=psfile] [, metafile=metafile] ; [, mindat=mindat] [, maxdat=maxdat] [, /linear] ; [, /sqroot] [, /logari] [, /invers] [, factor=factor] ; [, titles=titles] [, encaps=encaps] ; ; REQUIRED INPUTS: ; rawdat a struct containing the detector windows ; ; OPTIONAL INPUTS: ; ctable the colour table to be used for plotting the array ; ; KEYWORDS: ; psfile the name of the postscript file for a postscript plot ; metafile the name of the metafile file (currently not implemented!!) ; mindat lowest grayscale value for scaling of the plot ; maxdat highest grayscale value for scaling of the plot ; linear use linear scaling of the intensities ; sqroot use square root scaling of the intensities (default) ; linear use logarithmic scaling of the intensities ; invers invert the colour table (e.g. for PS plots) ; factor magnification factor for screen plotting (default is 3) ; titles string containing a title for the plot ; encaps make encapsulated PS (EPS) output ; ; OUTPUTS: ; a plot of the provided detector windows surrounded by multiple axes ; ; DESCRIPTION AND EXAMPLE: ; The programme has to be provided with a struct RAWDAT={data1:2d_array, ; data2:2d_array, ....} containing the detector frames from MIDI. Both ; dispersed data (e.g. photometry) or undispersed data (e.g. acquisition) ; is accepted. The programme then determines what kind of data is ; handled (dispersed: 'prism' or 'grism', undispersed: 'none') ; automatically from the pixel sizes of the detector windows. If MINDAT ; and MAXDAT are not provided, the programme determines these by ; rejecting the n brightest and m faintest pixel values. The two arrays ; are then rescaled and plotted. The lower x-axis and the left y-axis ; are always scaled with the pixel position on the detector, the ; alternative x- and y-axis are scaled with wavelength or angular offset ; (depending on the data type). For undispersed data ('none') the ; windows are plotted next to each other, for dispersed data ('prism' or ; 'grism') the windows are plotted on top of each other. The simplest ; way to run the programme is by calling for example: ; midi_pltwin, rawdat ; ; CALLED BY: ; none ; ; CALLING: ; midi_outdis ; midi_outfie ; ; MODIFICATION HISTORY: ; 2008-03-15 Written by Konrad R. Tristram ; 2008-05-14 Konrad R. Tristram: Change the way the reference is plotted. ; 2008-12-19 Konrad R. Tristram: Add keyword ENCAPS and change PS output. ; 2009-06-05 Konrad R. Tristram: Changed position of slit position to ; estimate derived from memo on 2008-12-22. ; 2009-06-30 Konrad R. Tristram: Added LLABEL and RLABEL to MIDI_OUTDIS. ; 2009-10-19 Konrad R. Tristram: Corrected error in the plot positioning. ; 2010-11-05 Konrad R. Tristram: Rewrite for arbitrary number of windows. ; 2010-12-06 Konrad R. Tristram: Added possibility to plot a title string. ; 2012-11-01 Konrad R. Tristram: Adjust PS font size for !p.font=0. ; PRO MIDI_PLTWIN, RAWDAT, CTABLE, PSFILE=PSFILE, METAFILE=METAFILE, $ MINDAT=MINDAT, MAXDAT=MAXDAT, LINEAR=LINEAR, SQROOT=SQROOT, $ LOGARI=LOGARI, INVERS=INVERS, FACTOR=FACTOR, TITLES=TITLES, $ ENCAPS=ENCAPS ; SET DEFAULT VALUES FOR SOME VARIABLES ;------------------------------------------------------------------------------- slitps = 1 ; show the slit position ; SAVE CURRENT PLOT SETTINGS AND TURN OFF MULTIPLE PLOTS ;------------------------------------------------------------------------------- p_old = !p d_old = !d.name !p.multi = 0 ; SET DEFAULT COLOR TABLE AND MAGNIFICATION (FOR SCREEN OUTPUT) ;------------------------------------------------------------------------------- if n_elements(ctable) eq 0 then ctable = 3 $ else ctable = byte(ctable[0]) < 45 if n_elements(factor) eq 0 then factor = 3 $ else factor = float(factor[0]) < 10 ; CATCH THE ATTEMPT TO MAKE A METAFILE ON A UNIX SYSTEM ;------------------------------------------------------------------------------- if !version.os_family ne 'Windows' and keyword_set(metafile) then begin print, 'Cannot produce METAFILES on Unix systems.' print, 'Producing PS output instead.' psfile = metafile endif ; FIND HOW MANY AND WHERE THE DATA WINDOWS ARE ;------------------------------------------------------------------------------- tagnam = tag_names(rawdat) winidx = where(strmatch(tagnam, 'DATA*')) if winidx[0] lt 0 then begin message, 'No detector data fields detected in the data. Exiting!', /cont return endif winidx = winidx[sort(tagnam[winidx])] winnum = n_elements(winidx) ; FIND OUT AT WHAT KIND OF DATA WE ARE LOOKING AT ;------------------------------------------------------------------------------- arrsiz = size(rawdat[0].(winidx[0])) case arrsiz[1] of 171: disper = 'prism' 261: disper = 'grism' 62: disper = 'none' else: disper = 'unknown' endcase ; FIND THE MINIMUM AND THE MAXIMUM OF THE DATA IF NOT PROVIDED ;------------------------------------------------------------------------------- alldat = float(rawdat[0].(winidx[0])) for i=1,winnum-1 do alldat = [[[alldat]], [[rawdat[0].(winidx[i])]]] srtdat = sort(alldat) if n_elements(mindat) ne 1 then mindat = alldat[srtdat[30]] if n_elements(maxdat) ne 1 then maxdat = alldat[srtdat[n_elements(alldat)-20]] ; NORMALISE THE DATA TO [0,1] AND ONLY SELECT THE FIRST FRAME ;------------------------------------------------------------------------------- alldat = (alldat - mindat) / (maxdat - mindat) > 0 < 1 ; RESCALE THE DATA AS SELECTED, DEFAULT IS SQUARE ROOT ;------------------------------------------------------------------------------- if keyword_set(linear) then begin alldat = alldat * 255 endif else if keyword_set(logari) then begin alldat = alog((alldat * exp(10)) + 1) * 25.5 endif else begin alldat = sqrt(alldat) * 255 endelse ; INVERT THE ARRAY TO BE PLOTTED IN ORDER TO GET BLACK ON WHITE PLOTS ;------------------------------------------------------------------------------- if keyword_set(invers) then alldat = 255 - alldat ; MAKE SETTINGS FOR THE SELECTED PLOTTING DEVICE ;------------------------------------------------------------------------------- if keyword_set(psfile) then begin ; SET CHARACTER SIZE, LOWER LEFT MARGIN OF PLOTS AND GAP BETWEEN PLOTS ;----------------------------------------------------------------------- if !p.font eq 0 then !p.charsize=0.65 else !p.charsize=0.75 posoff = [1000, 1000] posgap = 200 dscale = 1.0 ; DEFINE THE WITDH OF THE PLOT ;----------------------------------------------------------------------- x_size = 12.0 ; A&A intermediate plot x_size = 8.8 ; A&A single column plot x_size=12.0 ; DISTINGUISH BETWEEN DISPERSED AND NONDISPERSED DATA ;----------------------------------------------------------------------- if disper eq 'none' then begin ; NONDISPERSED WINDOWS ARE PLOTTED NEXT TO EACH OTHER ;--------------------------------------------------------------- factor = (1000 * x_size - 2*posoff[0] - posgap) / 2 / arrsiz[1] hlparr = [1, 0]#[1, 1] * (posgap + factor * arrsiz[1]) endif else begin ; DISPERSED WINDOWS ARE PLOTTED ON TOP OF EACH OTHER ;--------------------------------------------------------------- factor = (1000 * x_size - 2*posoff[0]) / arrsiz[1] hlparr = [0, 1]#[1, 1] * (posgap + factor * arrsiz[2]) endelse ; CALCULATE THE POSITIONS OF THE PLOTS ;----------------------------------------------------------------------- pltpos = [posoff, posoff + factor * arrsiz[1:2]] # replicate(1,winnum) for i=0,winnum-1 do pltpos[*,i] += i * hlparr y_size = pltpos[3, winnum-1]/1000 + 1.0 ; OPEN THE PS FILE AND SET THE PLOT AREA ;----------------------------------------------------------------------- set_plot, 'PS' psstri = 'systemdict /setdistillerparams known { ' + string(10b) + $ '<< /AutoFilterColorImages false /ColorImageFilter ' + $ string(10b) + '/FlateEncode >> setdistillerparams} if' device, file=psfile+'.ps', xsize=x_size, ysize=y_size, $ xoffset=10.4-0.5*x_size, yoffset=15.0-0.5*y_size, $ bits_per_pixel=8, /portrait, color=1, output=psstri, $ encaps=encaps endif else if keyword_set(metafile) then begin ; SET FONT, CHARACTER SIZE, LOWER LEFT MARGIN AND GAP BETWEEN PLOTS ;----------------------------------------------------------------------- !p.font=0 !p.charsize=0.80 posoff = [1000, 1000] posgap = 200 ; OPEN THE METAFILE FILE AND SET THE PLOT AREA ;----------------------------------------------------------------------- set_plot, 'METAFILE' device, file=metafile+'.emf', set_font='Arial*76', $ xsize=56.0, ysize=42.0 endif else begin ; SET CHARACTER SIZE, LOWER LEFT MARGIN OF PLOTS AND GAP BETWEEN PLOTS ;----------------------------------------------------------------------- !p.charsize=1.25 posoff = [50, 50] posgap = 20 dscale = 0.045 ; CALCULATE THE POSITIONS OF THE PLOTS ;----------------------------------------------------------------------- if disper eq 'none' then hlparr = [1, 0] else hlparr = [0, 1] pltpos = [posoff, posoff + factor * arrsiz[1:2]] # replicate(1,winnum) for i=0,winnum-1 do pltpos[*,i] += i * hlparr # [1, 1] * $ (posgap + factor * arrsiz[2]) ; OPEN THE DEVICE AND A NEW PLOT WINDOW ;----------------------------------------------------------------------- if !version.os_family eq 'Windows' then set_plot, 'WIN' $ else set_plot, 'X' window, xsize=pltpos[2, winnum-1]+50, ysize=pltpos[3,winnum-1]+50, /FREE device, decompose=0 endelse ; CALCULATE THE POSITIONS OF THE LABELS ;------------------------------------------------------------------------------- l1posi = [posoff[0]-dscale*700, (pltpos[3, winnum-1]+pltpos[1, 0])/2] l2posi = [pltpos[2, 0]+dscale*800, (pltpos[3, winnum-1]+pltpos[1, 0])/2] l3posi = [(pltpos[2, 0]+pltpos[0, winnum-1])/2, pltpos[3, winnum-1]+dscale*600] l4posi = [(pltpos[2, 0]+pltpos[0, winnum-1])/2, pltpos[3, winnum-1]+dscale*200] ; DISTINGUISH BETWEEN DISPERSED AND NONDISPERSED WINDOWS ;------------------------------------------------------------------------------- if disper eq 'none' then begin ; PLOT THE FIRST WINDOW ON THE LEFT AND LABEL IT ;----------------------------------------------------------------------- ; peakps = [[30.0, 32.0],[29.0, 33.0]] ; Position before 2006-01-01 ; peakps = [[29.0, 34.7],[28.0, 35.0]] ; Corrected position on 2008-04-16 ; peakps = [[29.0, 31.0],[28.0, 32.0]] ; Position before 2008-04-16 ; peakps = [[29.0, 31.7],[28.0, 32.0]] ; Position after 2008-04-16 peakps = [[30.0, 31.7],[29.0, 32.0]] ; Position after 2009-04-01 ; LOOP OVER THE DIFFERENT WINDOWS ;----------------------------------------------------------------------- for i=0,winnum-1 do begin ; PLOT THE CURRENT WINDOW AND LABEL IT ;--------------------------------------------------------------- if n_elements(titles) eq 0 then tlabel=1 if i eq 0 then llabel = 1 else llabel = 0 if i eq winnum-1 then rlabel = 1 else rlabel = 0 midi_outfie, alldat[*,*,i], pltpos[*,i], $ peakps[*,fix(2*i/winnum)], ctable=ctable, $ llabel=llabel, rlabel=rlabel, tlabel=tlabel xyouts, 5, 60, 'Window '+strtrim(i+1,2), charsiz=1.5*!p.charsize ; ADD A CROSS AT THE SOURCE POSITION ;--------------------------------------------------------------- loadct, 39, /silent plots, peakps[0,fix(2*i/winnum)], peakps[1,fix(2*i/winnum)], $ psym=7, color=100, symsize=1.5*!p.charsize ; SHOW THE POSITION OF THE 200 MICRON = 0.52 ARCSEC SLIT ;--------------------------------------------------------------- if keyword_set(slitps) then begin oplot, [1,1] # 27.5-fix(2*i/winnum), [-1,100], color=100 oplot, [1,1] # 32.5-fix(2*i/winnum), [-1,100], color=100 endif endfor endif else begin ; LOOP OVER THE DIFFERENT WINDOWS ;----------------------------------------------------------------------- for i=0,winnum-1 do begin ; PLOT THE CURRENT WINDOW AND LABEL IT ;--------------------------------------------------------------- if i eq 0 then blabel = 1 else blabel = 0 if i eq winnum-1 then tlabel = 1 else tlabel = 0 midi_outdis, alldat[*,*,winnum-i-1], pltpos[*,i], $ 15-fix(2*i/winnum), ctable=ctable, disper=disper, $ tlabel=tlabel, blabel=blabel, /llabel, /rlabel xyouts, 8, 31, 'Window '+strtrim(winnum-i,2), $ charsize=1.5*!p.charsize endfor ; PLOT THE AXES LABELS ;----------------------------------------------------------------------- xyouts, l1posi[0], l1posi[1], 'pixels on detector', $ align=0.5, orientation=90, /device xyouts, l2posi[0], l2posi[1], 'offset [arcsec]', $ align=0.5, orientation=90, /device if n_elements(titles) eq 0 then xyouts, l3posi[0], l3posi[1], $ 'wavelength !4k!X [!4l!Xm]', align=0.5, /device endelse ; PLOT A TITLE STRING IF ONE HAS BEEN PROVIDED ;------------------------------------------------------------------------------- if n_elements(titles) ne 0 then xyouts, l3posi[0], l3posi[1], titles, $ align=0.5, /device ; CLOSE THE FILES IF CREATING A PS FILE OR A METAFILE ;------------------------------------------------------------------------------- if keyword_set(psfile) or keyword_set(metafile) then device, /close ; RESTORE ORIGINAL PLOT SETTINGS ;------------------------------------------------------------------------------- set_plot, d_old !p = p_old END ;+ ; NAME: ; MIDI_PHOTOM ; ; PURPOSE: ; Do chopping for MIDI photometry data and extract the spectrum. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; result = midi_photom(phtdat [, mskfil=mskfil] [, sigma=sigma] ; [, /reset] [, /chop] [, /output] [, /nosky] ; [, porder=porder] [, phtbnd=phtbnd] ; [, skydis=skydis] [, skywid=skywid] ; [, skysmo=skysmo] [, psfile=psfile] ; [, /invers] [, titstr=titstr]) ; ; REQUIRED INPUTS: ; phtdat photometric data struct as produced by oirGetData ; e.g. phtdat = oirGetData('photA.fits') ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORDS: ; mskfil location of optional mask file, composed of directory ; and file name, e.g. ; mskfil = './minrtsMask_FIELD_PRISM_HIGH_SENS_SLIT.fits' ; sigma sigma for rejection of frames in guess of chopping state ; reset reset chopping flags to the ones determined by this program ; chop display variations in object flux due to chopping, high ; value is on target low value is on sky ; output display results in graphics window ; nosky don't do extra sky subtraction using the sky bands ; porder order of the polynomial which is fitted to the sky ; phtbnd two element with the y-range of the object band ; skydis distance of sky band from object band (default: 3 rows) ; skysmo width of boxcar to smooth sky (default: 10 columns) ; skywid width of sky band in pixels (default = 3) ; psfile the name of the postscript file for a postscript plot ; invers invert the colour table (e.g. for PS plots) ; titstr a title string to also plot ; ; OUTPUTS: ; result struct containg the sky-subtracted, compressed detector ; images and the spectra for both detector windows ; ; DESCRIPTION AND EXAMPLE: ; First some photometry data has to be loaded with oirGetData, for ; example phtdat = oirGetData('MIDI.2000-01-01T12:34:56.000.fits'). The ; simplest way to start the program is then by 'photometry, phtdat', ; which will take the flags in the imput struct to determine whether ; a frame is on sky or on target. All sky and all object frames are ; compressed and the sky subtracted from the object. An additional sky ; subtraction is done using two sky bands located SKYD away from the ; central part. In the default mode, the sky is simply the average of ; the two sky bands. When PORDER is set, a polynomial fit of the order ; PORDER is performed instead. The photometry of the object is ; determined from the masked area after sky subtraction. If no file ; name for a mask is provided, the mask is a top hat function, the ; position of which is given by PHTBND. The result is stored in a ; struct returned by the function. Options include to visualise the ; difference between object and sky frames and a reset of the flags in ; the input file from this information. The result can be displayed ; with the option /OUTPUT. ; ; CALLED BY: ; midi_phtfil ; ; CALLING: ; midi_outdis ; midi_wavaxi ; ; MODIFICATION HISTORY: ; 2005-09-02 Written by Konrad R. Tristram ; 2006-03-28 Konrad R. Tristram: Change in documentation, setting of the ; photometry area and yrange of the plot. ; 2006-10-10 Konrad R. Tristram: Added possibility to plot into PS file. ; 2007-12-10 Konrad R. Tristram: Added parameter SKYSMO. ; 2008-03-15 Konrad R. Tristram: Completely rewritten, multiple changes. ; 2008-05-29 Konrad R. Tristram: Mask & sky now 1 pixel higher in win 2. ; 2008-08-08 Konrad R. Tristram: Added polynomial fit for sky background. ; 2008-11-06 Konrad R. Tristram: Added polynomial fit for sky background. ; 2008-12-11 Konrad R. Tristram: Inhibit compression for PS output. ; 2008-12-11 Konrad R. Tristram: Don't normalise counts by mask. ; 2009-06-30 Konrad R. Tristram: Added LLABEL and RLABEL to MIDI_OUTDIS. ; 2009-06-30 Konrad R. Tristram: Fixed error in smoothing part. ; 2012-04-03 Konrad R. Tristram: Changed counters in loops to LONG. ; 2012-04-03 Konrad R. Tristram: Renamed keyword TITSTR to TITLES. ; 2012-11-01 Konrad R. Tristram: Adjust PS font size for !p.font=0. ; FUNCTION MIDI_PHOTOM, PHTDAT, MSKFIL=MSKFIL, SIGMA=SIGMA, RESET=RESET, $ CHOP=CHOP, OUTPUT=OUTPUT, NOSKY=NOSKY, PORDER=PORDER, $ PHTBND=PHTBND, SKYDIS=SKYDIS, SKYWID=SKYWID, $ SKYSMO=SKYSMO, PSFILE=PSFILE, INVERS=INVERS, $ TITLES=TITLES ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- if (n_elements(mskfil) ne 1) then mskfil = '' if (n_elements(sigma) ne 1) then sigma = 3.0 if (n_elements(skydis) ne 1) then skydis = 4 else skydis = fix(skydis)+1 if (n_elements(skywid) ne 1) then skyw = 2 else skyw = skywid - 1 if (n_elements(skysmo) ne 1) then skysmo = 10 else skysmo = byte(skysmo) if (n_elements(phtbnd) ne 2) then begin ymin = 12b ; minimum line of photometry area ymax = 21b ; maximum line of photometry area endif else begin ymin = byte(phtbnd[0]) < 40 ymax = byte(phtbnd[1]) < 40 endelse xmin = 40 ; minimum column of photometry area xmax = 130 ; maximum column of photometry area ; SAVE THE CURRENT PLOT SETTINGS AND TURN OFF MULTIPLE PLOTS ;------------------------------------------------------------------------------- p_old = !p d_old = !d.name p_multi = !p.multi !p.multi = 0 ; FIND OUT IF WE ARE REALLY LOOKING AT PRISM DATA ;------------------------------------------------------------------------------- arrsiz = size(phtdat[0].data1) if ((arrsiz[1] ne 151) and (arrsiz[1] ne 171)) or (arrsiz[2] ne 41) then begin print, 'Only know how to deal with PRISM data. Returning.' return, -1 endif else disper = 'prism' ; GO INTO CHOPPING STATUS DETERMINATION PART ;------------------------------------------------------------------------------- if keyword_set(CHOP) or keyword_set(RESET) then begin ; MAKE A SIMPLE SKY REMOVAL AND ENHANCE THE CONTRAST BY THE TAKING DIFFERENCE ;------------------------------------------------------------------------------- data1 = phtdat.data1-smooth(phtdat.data1, [0,0,40]) data2 = phtdat.data2-smooth(phtdat.data2, [0,0,40]) data = data1 + data2 ; FOR EVERY FRAME DETERMINE OBJECT FLUX (OBJECT-SKY) AND NOISE IN SKY BAND ;------------------------------------------------------------------------------- difference = fltarr(n_elements(data[1,1,*])) variations = fltarr(n_elements(data[1,1,*])) for i=0L,n_elements(data[1,1,*])-1 do begin obj = 1.0*data[xmin:xmax,ymin:ymax,i] sky = [(data[xmin:xmax,ymin-skydis-skyw:ymin-skydis,i]), $ data[xmin:xmax,ymin+skydis:ymin+skydis+skyw,i]] difference[i]=total(obj)/n_elements(obj)-total(sky)/n_elements(sky) variations[i]=sigma*stdev(sky)/sqrt(n_elements(sky)) endfor ; SMOOTH OBJECT FLUX AND NOISE ESTIMATE BECAUSE THEY ARE VERY "SPIKY" ;------------------------------------------------------------------------------- difference = smooth(difference, 6) variations = smooth(variations, 6) ; PLOT VARIATION OF OBJECT FLUX DUE TO CHOPPING ;------------------------------------------------------------------------------- if keyword_set(CHOP) or keyword_set(OUTPUT) then begin if !version.os_family eq 'Windows' then set_plot, 'WIN' $ else set_plot, 'X' device, decompose=0 window, /FREE, xsize=1400, ysize=400 plot, difference, /NODATA loadct, 33, /silent oplot, difference, color=162 oplot, variations, color=230 oplot, -variations, color=230 xyouts, 25, 200, 'Check Chopping', charsize=2, alignment=0.5, $ ORIENTATION=90, /DEVICE ; window, /FREE & plot, fft(difference) endif endif ; OVERWRITE CHOPPING STATUS FLAGS IN INPUT STRUCT WITH NEW CHOPPING STATUS ;------------------------------------------------------------------------------- if keyword_set(RESET) then begin for i=0L,n_elements(data[1,1,*])-1 do begin if difference[i] lt variations[i] then targtype = 'S ' $ else if difference[i] gt variations[i] then targtype = 'T ' $ else targtype = 'U ' phtdat[i].tartyp1=targtype phtdat[i].tartyp2=targtype endfor endif ; FIND SKY AND OBJECT FRAMES AND SUBTRACT SKY FROM OBJECT FOR BOTH WINDOWS ;------------------------------------------------------------------------------- obj = where(phtdat.tartyp1 eq "T ") sky = where(phtdat.tartyp1 eq "S ") result1 = float(total(phtdat[obj].data1, 3))/n_elements(obj) - $ float(total(phtdat[sky].data1, 3))/n_elements(sky) obj = where(phtdat.tartyp2 eq "T ") sky = where(phtdat.tartyp2 eq "S ") result2 = float(total(phtdat[obj].data2, 3))/n_elements(obj) - $ float(total(phtdat[sky].data2, 3))/n_elements(sky) ; CHECK FOR MASK AND GET MASK DATA OR CREATE OWN TOPHAT MASK ;------------------------------------------------------------------------------- info=file_info(mskfil) if info.regular then begin print, "Mask provided was found!" mask=oirgetdata(mskfil) endif else begin print, "No mask was found, using tophat mask." mask={data1:make_array(171, 41, /FLOAT, value=0.0), $ data2:make_array(171, 41, /FLOAT, value=0.0)} mask.data1[*,ymin:ymax] = 1.0 mask.data2[*,ymin+1:ymax+1] = 1.0 endelse ; DO EXTRA SKY SUBTRACTION BY A POLYNOMIAL FIT OF COLUMNS IN Y DIRECTION ;------------------------------------------------------------------------------- if n_elements(porder) eq 1 then begin porder = porder > 0 < 10 sky = result1 ; DO THE SUBTRACTION FOR WINDOW 1 ;----------------------------------------------------------------------- sel = [ymin-skydis-indgen(skyw+1), ymax+skydis+indgen(skyw+1)] yps = indgen(arrsiz[2]) for i=0L,arrsiz[1]-1 do begin $ col = result1[(i-skysmo/2>0):(i+skysmo/2<(arrsiz[1]-1)),*] col = total(col,1) / n_elements(col[*,0]) sky[i,*] = poly(yps, poly_fit(sel, col[sel], porder)) endfor result1 = result1 - sky ; DO THE SUBTRACTION FOR WINDOW 2 WHICH IS ONE PIXEL HIGHER ;----------------------------------------------------------------------- sel = sel; + 1 for i=0L,arrsiz[1]-1 do begin $ col = result2[(i-skysmo/2>0):(i+skysmo/2<(arrsiz[1]-1)),*] col = total(col,1) / n_elements(col[*,0]) sky[i,*] = poly(yps, poly_fit(sel, col[sel], porder)) endfor result2 = result2 - sky ; OTHERWISE DO FAST LINEAR FIT BY SIMPLY AVERAGING THE TOP AND BOTTOM SKY BANDS ;------------------------------------------------------------------------------- endif else if not keyword_set(NOSKY) then begin yps = replicate(1.0, n_elements(result1[0,*])) ; DO THE SUBTRACTION FOR WINDOW 1 ;----------------------------------------------------------------------- sky = 0.5*(result1[*,ymin-skydis-skyw:ymin-skydis] + $ result1[*,ymax+skydis:ymax+skydis+skyw]) sky = smooth(total(sky,2)/n_elements(sky[1,*]), skysmo) result1 = result1 - sky # yps ; DO THE SUBTRACTION FOR WINDOW 2 WHICH IS ONE PIXEL HIGHER ;----------------------------------------------------------------------- sky = 0.5*(result2[*,ymin-skydis-skyw+1:ymin-skydis+1] + $ result2[*,ymax+skydis+1:ymax+skydis+skyw+1]) sky = smooth(total(sky,2)/n_elements(sky[1,*]), skysmo) result2 = result2 - sky # yps endif ; APPLY MASK TO DATA AND DETERMINE THE SPECTRUM ;------------------------------------------------------------------------------- spectrum1 = total(result1*mask.data1, 2);/total(mask.data1, 2) spectrum2 = total(result2*mask.data2, 2);/total(mask.data2, 2) ; START OUTPUT IF KEYWORD 'OUTPUT' IS SET ;------------------------------------------------------------------------------- if keyword_set(OUTPUT) then begin ; MAKE SETTINGS FOR THE SELECTED PLOTTING DEVICE ;----------------------------------------------------------------------- if keyword_set(psfile) then begin ; SET CHARACTER SIZE, LOWER LEFT MARGIN OF PLOTS AND GAP ;--------------------------------------------------------------- if !p.font eq 0 then !p.charsize=0.65 else !p.charsize=0.75 posoff = [1000, 1000] posgap = 200 ; DEFINE THE WITDH OF THE PLOT ;--------------------------------------------------------------- x_size = 12.0 ; A&A intermediate plot ;x_size = 8.8 ; A&A single column plot ; CALCULATE THE POSITIONS OF THE PLOTS AND THE LABELS ;--------------------------------------------------------------- factor = (1000 * x_size - 2*posoff[0]) / arrsiz[1] hlparr = [0, 1] posit2 = [posoff, posoff[0] + factor * arrsiz[1], $ posoff[1] + factor * arrsiz[2]] posit1 = posit2 + [0, 1]#[1, 1] * (posgap + factor * arrsiz[2]) posit3 = [posoff[0]-700, (posit2[3]+posit1[1])/2] posit4 = [posit1[2]+800, (posit2[3]+posit1[1])/2] posit6 = posit1 posit6[1] = posit1[3] + 1000 posit6[3] = posit6[1] + (posit6[2]-posit6[0]) / 2 posit5 = [(posit6[2]+posit6[0])/2, posit6[3]+600] y_size = posit6[3]/1000 + 1.0 ; CALCULATE THE POSITIONS FOR THE TITLE STRING ;--------------------------------------------------------------- if n_elements(titstr) eq 1 then begin posit7 = [200, posit5[1] + 500] y_size = y_size + 0.5 endif ; OPEN THE PS FILE AND SET THE PLOT AREA ;--------------------------------------------------------------- set_plot, 'PS' psstri = 'systemdict /setdistillerparams known { ' + $ string(10b) + '<< /AutoFilterColorImages false ' + $ '/ColorImageFilter ' + string(10b) + $ '/FlateEncode >> setdistillerparams} if' device, file=psfile+'.ps', xsize=x_size, ysize=y_size, $ xoffset=10.4-0.5*x_size, yoffset=15.0-0.5*y_size, $ bits_per_pixel=8, /portrait, color=1, output=psstri endif else if keyword_set(metafile) then begin ; OPEN METAFILE AND SET PLOT AREA ;--------------------------------------------------------------- set_plot, 'METAFILE' device, file=metafile+'.emf', set_font='Arial*76', $ xsize=56.0, ysize=42.0 ; SET FONT PROPERTIES AND POSITIONS FOR OUTPUT ;--------------------------------------------------------------- !p.font=0 !p.charsize=0.80 endif else begin ; SET CHARACTER SIZE, LOWER LEFT MARGIN, GAP AND MAGNIFICATION ;--------------------------------------------------------------- !p.charsize=1.25 posoff = [50, 40] posgap = 10 factor = 3 ; CALCULATE THE POSITIONS OF THE PLOTS AND OF THE LABELS ;--------------------------------------------------------------- if disper eq 'none' then hlparr = [1, 0] else hlparr = [0, 1] posit2 = [posoff, posoff + factor * arrsiz[1:2]] posit1 = posit2 + hlparr # [1, 1] * (10 + factor * arrsiz[2]) posit3 = [posoff[0]-30, (posit2[3]+posit1[1])/2] posit4 = [posit1[2]+35, (posit2[3]+posit1[1])/2] posit6 = posit1 posit6[1] = posit1[3] + 50 posit6[3] = posit6[1] + (posit6[2]-posit6[0]) / 2 posit5 = [(posit6[2]+posit6[0])/2, posit6[3]+30] posit7 = [10, posit5[1]] y_size = posit6[3] + 50 ; CALCULATE THE POSITIONS FOR THE TITLE STRING ;--------------------------------------------------------------- if n_elements(titstr) eq 1 then begin posit7 = [10, posit5[1] + 20] y_size = y_size + 20 endif ; OPEN THE DEVICE AND A NEW PLOT WINDOW ;--------------------------------------------------------------- if !version.os_family eq 'Windows' then set_plot, 'WIN' $ else set_plot, 'X' window, xsize=posit1[2]+40, ysize=y_size, /FREE device, decompose=0 endelse ; SET THE COLOURS AN THE COLOUR TABLE ;----------------------------------------------------------------------- colour1 = 200 colour2 = 230 ctable = 3 ; PLOT THE SPECTRA AND THE SKY IN THE UPPER PART ;----------------------------------------------------------------------- loadct, 39, /silent plot, spectrum1, xrange=[-0.5, arrsiz[1]-0.5], xstyle=9, $ yrange=[min([spectrum1,spectrum2]),max([spectrum1,spectrum2])], $ ystyle=0, ytitle="counts in ADU", $ position=posit6, /NODATA, /DEVICE, /NOERASE oplot, [-1000,1000], [0,0], linestyle=1 oplot, smooth(total(result1[*,ymin-skydis-skyw:ymin-skydis] + $ result1[*,ymax+skydis:ymax+skydis+skyw], 2) / $ (2*skyw+2), skysmo), color=colour1, linestyle=1 oplot, smooth(total(result2[*,ymin-skydis-skyw+1:ymin-skydis+1] + $ result2[*,ymax+skydis+1:ymax+skydis+skyw+1], 2) / $ (2*skyw+2), skysmo), color=colour2, linestyle=1 oplot, spectrum1, color=colour1 oplot, spectrum2, color=colour2 ; ADD THE WAVELENGTH AXES ON TOP OF THE PLOT ;----------------------------------------------------------------------- midi_wavaxi, disper, /pltlab if !p.font gt -1 then strvar = ['!9l!X [!9m!Xm]'] $ else strvar = ['!4k!X [!4l!Xm]'] xyouts, posit5[0], posit5[1], 'wavelength '+strvar, align=0.5, /device ; PLOT THE TITLE STRING IF PROVIDED ;----------------------------------------------------------------------- if n_elements(titstr) eq 1 then begin xyouts, posit7[0], posit7[1], titstr, align=0.0, $ charsize=1.25*!p.charsize, /device endif ; FIND THE MINIMUM AND THE MAXIMUM OF THE DATA ;----------------------------------------------------------------------- if n_elements(porder) eq 1 then begin alldat = [result1[*,ymin-skydis-skyw+0:ymax+skydis+skyw+0], $ result2[*,ymin-skydis-skyw+1:ymax+skydis+skyw+1]] alldat = float(alldat) endif else alldat = [result1, result2] alldat = float(alldat) srtdat = sort(alldat) mindat = alldat[srtdat[40]] maxdat = alldat[srtdat[n_elements(alldat)-30]] ; NORMALISE DATA TO [0,1] ;----------------------------------------------------------------------- data1 = (float(result1) - mindat) / (maxdat - mindat) > 0 < 1 data2 = (float(result2) - mindat) / (maxdat - mindat) > 0 < 1 ; RESCALE THE DATA LINEARLY TO 1 BYTE COLOUR DEPTH ;----------------------------------------------------------------------- data1 = sqrt(data1) * 255 data2 = sqrt(data2) * 255 ; TO SHOW THE TRUE SKY BANDS UNCOMMENT THE FOLLOWING LINES ;----------------------------------------------------------------------- ;data1[*,ymin-skydis-skyw:ymin-skydis] = 255 ;data1[*,ymax+skydis:ymax+skydis+skyw] = 255 ;data2[*,ymin-skydis-skyw+1:ymin-skydis+1] = 255 ;data2[*,ymax+skydis+1:ymax+skydis+skyw+1] = 255 ; TO SHOW THE TRUE OBJECT BAND UNCOMMENT THE FOLLOWING LINES ;----------------------------------------------------------------------- ;data1[*,ymin:ymax] = 255 ;data2[*,ymin+1:ymax+1] = 255 ; INVERT THE ARRAY TO BE PLOTTED IN ORDER TO GET BLACK ON WHITE PLOTS ;----------------------------------------------------------------------- if keyword_set(invers) then begin data1 = 255 - data1 data2 = 255 - data2 colour3 = 100 endif else colour3 = 100 ; PLOT OUT THE TWO WINDOWS ON TOP OF EACH OTHER AND LABEL THEM ;----------------------------------------------------------------------- midi_outdis, data1, posit1, 14, ctable=ctable, disper=disper, $ /tlabel, /llabel, /rlabel loadct, 39, /silent oplot, [140, 145], [ymin, ymin]-0.5, color=colour3, thick=2 oplot, [140, 145], [ymin-skydis, ymin-skydis]+0.5, color=colour3 oplot, [140, 145], [ymin-skydis-skyw, ymin-skydis-skyw]-0.5, $ color=colour3 oplot, [140, 145], [ymax, ymax]+0.5, color=colour3, thick=2 oplot, [140, 145], [ymax+skydis, ymax+skydis]-0.5, color=colour3 oplot, [140, 145], [ymax+skydis+skyw, ymax+skydis+skyw]+0.5, $ color=colour3 xyouts, 8, 31, 'Window 1', charsize=1.5*!p.charsize, color=colour1 midi_outdis, data2, posit2, 15, ctable=ctable, disper=disper, $ /blabel, /llabel, /rlabel loadct, 39, /silent oplot, [140, 145], [ymin+1, ymin+1]-0.5, color=colour3, thick=2 oplot, [140, 145], [ymin-skydis+1, ymin-skydis+1]+0.5, color=colour3 oplot, [140, 145], [ymin-skydis-skyw+1, ymin-skydis-skyw+1]-0.5, $ color=colour3 oplot, [140, 145], [ymax+1, ymax+1]+0.5, color=colour3, thick=2 oplot, [140, 145], [ymax+skydis+1, ymax+skydis+1]-0.5, color=colour3 oplot, [140, 145], [ymax+skydis+skyw+1, ymax+skydis+skyw+1]+0.5, $ color=colour3 xyouts, 8, 31, 'Window 2', charsize=1.5*!p.charsize, color=colour2 ; PLOT THE AXES LABELS ;----------------------------------------------------------------------- xyouts, posit3[0], posit3[1], 'pixels on detector', $ align=0.5, orie=90, /device xyouts, posit4[0], posit4[1], 'offset [arcsec]', $ align=0.5, orie=90, /device ; CLOSE THE FILES IF CREATING A PS FILE OR A METAFILE ;----------------------------------------------------------------------- if keyword_set(psfile) or keyword_set(metafile) then device, /close endif ; RESTORE THE ORIGINAL PLOT SETTINGS ;------------------------------------------------------------------------------- set_plot, d_old !p = p_old !p.multi = p_multi ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, {data1:result1, data2:result2, spectrum1:spectrum1, spectrum2:spectrum2} END ;+ ; NAME: ; MIDI_PHTFIL ; ; PURPOSE: ; Reduce a MIDI photometry data file. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; result = midi_phtfil(phtfil [, mskfil=mskfil] [, /nosky] ; [, phtbnd=phtbnd] [, skydis=skydis] ; [, psfile=psfile] [, /invers]) ; ; REQUIRED INPUTS: ; phtfil name of the photometry file ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORDS: ; mskfil location of optional mask file, composed of directory ; and file name, e.g. ; mskfil = './minrtsMask_FIELD_PRISM_HIGH_SENS_SLIT.fits' ; nosky don't do extra sky subtraction using the sky bands ; phtbnd two element with the y-range of the object band ; skydis distance of sky band from object band (default: 3 rows) ; psfile the name of the postscript file for a postscript plot ; invers invert the colour table (e.g. for PS plots) ; ; OUTPUTS: ; result struct containg the sky-subtracted, compressed detector ; images and the spectra for both detector windows ; ; DESCRIPTION AND EXAMPLE: ; This function is a wrapper function for loading the photometric data ; and then calling the reduction and plotting function MIDI_PHOTOM. For ; a detailed description of the options and the output see there. The ; simplest way to call the function is for example: ; result=midi_phtfil('./MIDI.2000-01-01T12:34:56.000.fits') ; ; CALLED BY: ; none ; ; CALLING: ; midi_photom ; ; MODIFICATION HISTORY: ; 2008-03-15 Written by Konrad R. Tristram ; 2012-04-03 Konrad R. Tristram: Renamed variable TITSTR to TITLES. ; FUNCTION MIDI_PHTFIL, PHTFIL, MSKFIL=MSKFIL, NOSKY=NOSKY, PHTBND=PHTBND, $ SKYDIS=SKYDIS, PSFILE=PSFILE, INVERS=INVERS ; CHECK IF THE FILE EXISTS ;------------------------------------------------------------------------------- if not total(file_test(strsplit(phtfil, ' ', /extract))) then begin print, 'File provided does not exist. Returning.' return, -1 endif ; GET THE SHUTTER POSITION AND CHECK IF DEALING WITH PHOTOMETRY DATA ;------------------------------------------------------------------------------- shuttr = oirGetKey(phtfil[0], 'INS SHUT NAME') if shuttr[0] eq 'AOPEN ' then shuttr = 'A ' else $ if shuttr[0] eq 'BOPEN ' then shuttr = 'B ' else begin print, 'No photometry data or something wrong with the file. Returning.' return, -1 endelse ; LOAD THE DATA ;------------------------------------------------------------------------------- phtdat = oirgetdata(phtfil) print, 'Data successfully loaded.' ; GET THE OBJECT NAME ;------------------------------------------------------------------------------- objnam = strtrim(oirGetKey(phtfil[0], 'OBS TARG NAME'), 2) titstr = 'Photometry ' + shuttr + 'of ' + objnam[0] ; RUN MIDI_PHOTOM ;------------------------------------------------------------------------------- result = midi_photom(phtdat, mskfil=mskfil, nosky=nosky, phtbnd=phtbnd, $ skydis=skydis, psfile=psfile, invers=invers, $ titles=titles, /output) ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, result END ;+ ; NAME: ; MIDI_PLTFRI ; ; PURPOSE: ; Show a video of the delay line position and the fringe on the detector. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; midi_pltfri, ifdata [, frames] [, delay] [, msknam=msknam] ; [, tsmoot=tsmoot] [, /davera] [, factor=factor] ; ; REQUIRED INPUTS: ; ifdata interferometric data struct as produced by oirGetData ; e.g. ifdata = oiGetData('MIDI.2000-01-01T12:34:56.000.fits') ; ; OPTIONAL INPUTS: ; frames start and end of the video in frame numbers ; delay delay between the display of every frame (in seconds) ; ; KEYWORDS: ; msknam name of a mask file that will be applied to the data ; tsmoot width of the boxcar for high pass filtering of the data ; davera additionally subtract the median from every frame ; factor magnification factor for screen plotting (default is 3) ; ; OUTPUTS: ; None. ; ; DESCRIPTION AND EXAMPLE: ; First some fringe tracking data has to be loaded with oirGetData, for ; example ifdata = oirGetData('MIDI.2000-01-01T12:34:56.000.fits'). The ; simplest way to start the program is then by 'fringedisplay, ifdata', ; which will then display the optical path distances (both the delay line ; and the piezos) in the upper half of the window and the detector data ; in the lower part of the window. The sum and the difference of the two ; windows are also plotted. ; Depending on the value of the keyword DELAY the programme runs in ; interactive mode (DELAY is a string variable) or in automatic mode ; (DELAY is undefined or a number). When a number is typed in the ; interactive mode, the programme switches to automatic mode using the ; number as the delay, if 'f' is typed, the programme finishes, by ; setting DELAY = 0.0. ; ; CALLED BY: ; none ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2005-09-02 Written by Konrad R. Tristram ; 2009-06-30 Konrad R. Tristram: Completely rewritten, multiple changes. ; 2012-04-03 Konrad R. Tristram: Changed counters in loops to LONG. ; PRO MIDI_PLTFRI, IFDATA, FRAMES, DELAY, MSKNAM=MSKNAM, TSMOOT=TSMOOT, $ DAVERA=DAVERA, FACTOR=FACTOR ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- method = 0 numpar = n_params() if (numpar lt 2) then frames = [1000, 2000] if not total(size(delay, /type) eq [1,2,3,4,5,7]) then delay = 0.005 $ else delay = delay[0] if n_elements(tsmoot) lt 1 then tsmoot = 20 else tsmoot = tsmoot > 0 < 500 if n_elements(msknam) ne 1 then msknam='' if n_elements(factor) eq 0 then factor = 3 $ else factor = float(factor[0]) < 10 ; SET CHARACTER SIZE ;------------------------------------------------------------------------------- chrsiz = !p.charsize !p.charsize = 1.25 ; CORRECT THE FRAME SELECTION SO THAT PROGRAMME WILL NOT CRASH ;------------------------------------------------------------------------------- frames[0] = frames[0] > tsmoot < (n_elements(ifdata) - tsmoot - 2) frames[1] = frames[1] > (frames[0] + 1) < (n_elements(ifdata) - tsmoot - 1) ; GET FRINGE DATA ;------------------------------------------------------------------------------- fringA = float(ifdata[(frames[0]-tsmoot):(frames[1]+tsmoot)].data1) fringB = float(ifdata[(frames[0]-tsmoot):(frames[1]+tsmoot)].data2) numele = n_elements(fringA[0,0,*]) xdimen = n_elements(fringA[*,0,0]) ydimen = n_elements(fringA[0,*,0]) ; CHECK FOR MASK AND GET MASK DATA OR CREATE TRANSPARENT MASK ;------------------------------------------------------------------------------- tmpvar = file_info(msknam) if tmpvar.exists then begin print, "Mask provided was found!" tmpvar = oirgetdata(msknam) maskwA = float(tmpvar.data1) maskwB = float(tmpvar.data2) endif else begin print, "No mask was found!" maskwA = make_array(xdimen, ydimen, /float, value=1.0) maskwB = make_array(xdimen, ydimen, /float, value=1.0) endelse ; USE TWO DIFFERENT METHODS TO REMOVE THE BACKGROUND ;------------------------------------------------------------------------------- if method then begin ; DO SIMPLE SKY SUBTRACTION (SUBTRACT OVERALL AVERAGE) AND MASKING ;----------------------------------------------------------------------- tmpvar = total(fringA,3) / numele for i=0L,numele-1 do fringA[*,*,i] = (fringA[*,*,i]-tmpvar) * maskwA tmpvar = total(fringB,3) / numele for i=0L,numele-1 do fringB[*,*,i] = (fringB[*,*,i]-tmpvar) * maskwB endif else begin ; DO HIGH PASS FILTERING IN TEMPORAL DIRECTION AND MASKING ;----------------------------------------------------------------------- fringA = fringA - smooth(fringA, [0,0,tsmoot]) for i=0L,numele-1 do fringA[*,*,i] = fringA[*,*,i] * maskwA fringB = fringB - smooth(fringB, [0,0,tsmoot]) for i=0L,numele-1 do fringB[*,*,i] = fringB[*,*,i] * maskwB endelse ; OPTIONALLY SUBTRACT THE MEDIAN FROM EACH FRAME ;------------------------------------------------------------------------------- if keyword_set(davera) then begin for i=0L,numele-1 do begin fringA[*,*,i] = fringA[*,*,i] - median(fringA[*,*,i]) fringB[*,*,i] = fringB[*,*,i] - median(fringB[*,*,i]) endfor endif ; CALCULATE THE SUM AND THE DIFFERENCE OF THE TWO WINDOWS ;------------------------------------------------------------------------------- fring1 = fringA + shift(fringB, 0, -1, 0) fring2 = fringA - shift(fringB, 0, -1, 0) ; MAKE RESULTING DATA SUITABLE FOR DISPLAYING, I.E. SET CONTRAST ;------------------------------------------------------------------------------- tmpvar = 7 * stddev(fring2) fringA=((fringA + tmpvar) * 127 / tmpvar) > 0 < 255 fringB=((fringB + tmpvar) * 127 / tmpvar) > 0 < 255 fring1=((fring1 + tmpvar) * 127 / tmpvar) > 0 < 255 fring2=((fring2 + tmpvar) * 127 / tmpvar) > 0 < 255 ; CALCULATE THE POSITIONS OF THE PLOTS AND OF THE LABELS ;------------------------------------------------------------------------------- posoff = [40, 40] posgap = 10 posit2 = [posoff, posoff + factor * [xdimen, ydimen]] posit1 = posit2 + [0, 1] # [1, 1] * (posgap + factor * ydimen) posit3 = posit1 + [1, 0] # [1, 1] * (posgap + factor * xdimen) posit4 = posit2 + [1, 0] # [1, 1] * (posgap + factor * xdimen) posit5 = [posoff[0]-30, (posit2[3]+posit1[1])/2] posit6 = [posit3[2]+35, (posit2[3]+posit1[1])/2] posit7 = [posoff[0], posit3[3]+ 40, posit3[2], 4*factor*ydimen+90] ; OPEN NEW WINDOW AND DRAW GRAPH ;------------------------------------------------------------------------------- window, /FREE, xsize=2*factor*xdimen+90, ysize=4*factor*ydimen+100 plot, [0,1], xrange=[frames[0],frames[1]], yrange=[-7.0,3.0], ystyle=1, $ xticklen=0.04, yticklen=0.01, xtitle='frame number', $ ytitle='OPD [!4l!Xm]', position=posit7, /device, /nodata ; PLOT THE LABELS FOR THE Y-AXES ;------------------------------------------------------------------------------- xyouts, posit5[0], posit5[1], 'pixels on detector', align=0.5, orie=90, /device xyouts, posit6[0], posit6[1], 'offset [arcsec]', align=0.5, orie=90, /device ; LOOP OVER THE FRAMES TO BE PLOTTED ;------------------------------------------------------------------------------- for i=0,frames[1]-frames[0] do begin ; PLOT THE DETECTOR DATA INCLUDING THE CORRESPONDING LABELS ;----------------------------------------------------------------------- midi_outdis, fringA[*,*,tsmoot+i], posit1, 14, ctable=0, $ disper='prism', /llabel loadct, 39, /silent xyouts, 6, 3, 'A', charsize=2, color=220, charthick=2 midi_outdis, fringB[*,*,tsmoot+i], posit2, 14, ctable=0, $ disper='prism', /blabel, /llabel loadct, 39, /silent xyouts, 6, 3, 'B', charsize=2, color=220, charthick=2 midi_outdis, fring1[*,*,tsmoot+i], posit3, 14, ctable=0, $ disper='prism', /rlabel loadct, 39, /silent xyouts, 6, 3, 'A + B', charsize=2, color=220, charthick=2 midi_outdis, fring2[*,*,tsmoot+i], posit4, 14, ctable=0, $ disper='prism', /blabel, /rlabel loadct, 39, /silent xyouts, 6, 3, 'A - B', charsize=2, color=220, charthick=2 ; RECREATE THE COORDINATE SYSTEM FOR THE OPD PLOT ;----------------------------------------------------------------------- loadct, 39, /silent plot, [0,1], xrange=[frames[0],frames[1]], yrange=[-7.0,3.0], $ xstyle=12, ystyle=13, position=posit7, /device, /noerase ; PLOT THE CURRENT OPD ;----------------------------------------------------------------------- oplot, [frames[0]+i], 1000.*[-20*ifdata[frames[0]+i].localopd[0]], $ psym=3, color=90 oplot, [frames[0]+i], 1000.*[ifdata[frames[0]+i].opd[1]], $ psym=3, color=250 ; WAIT FOR USER INPUT OR SOME PREDEFINED PERIOD BEFORE THE NEXT FRAME ;----------------------------------------------------------------------- if size(delay, /type) eq 7 then begin read, 'Pess ENTER to continue, "f" to finish!', delay if float(delay) ne 0.0 then delay = float(delay) $ else if delay eq 'f' then delay = 0.0 endif else wait, delay endfor ; RESET CHARACTER SIZE ;------------------------------------------------------------------------------- !p.charsize = chrsiz END ;+ ; NAME: ; MIDI_GAUDIS ; ; PURPOSE: ; Fit a Gaussian to every channel in a dispersed MIDI window. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; result = midi_gaudis(detwin [, rrange=rrange] [, smtval=smtval] ; [, /output] [, /quiet] ; ; REQUIRED INPUTS: ; detwin 2D array corresponding to a MIDI detector window ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORDS: ; rrange the range of rows used for the fitting ; smtval two element array with the smoothing values in y and y ; output plot a comparison between the data and the fit ; invers suppress all informational messages ; ; OUTPUTS: ; result struct containg the resulting fit parameters and the chi^2 ; ; DESCRIPTION AND EXAMPLE: ; This function fits a 1D Gaussian distribution to every column in a ; dispersed MIDI window (e.g. to photometries for PRISM or for GRISM ; data) using the function GAUSSFIT implemented in IDL. For the fit ; only the rows specified by the parameter RRANGE are used. If the ; keyword /OUTPUT is set then a comparison between the data and the ; fits will be also plotted. The resulting parameters as a function of ; column (the scaling factor for the Gaussian, the position of the peak ; and the FWHM) and their errors as well as the chi^2 are returned in ; an array of struct. The function is called by: ; result = midi_gaudis(detwin) ; ; CALLED BY: ; midi_mskfit ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2008-10-13 Written by Konrad R. Tristram ; 2008-11-19 K. Tristram: Added the possibility to smooth DETWIN. ; 2009-04-26 K. Tristram: Guess the row where the spectrum is located. ; 2010-11-02 K. Tristram: Changed array of struct to struct of arrays. ; 2011-07-13 K. Tristram: Adapted to IDL 8.1. ; FUNCTION MIDI_GAUDIS, DETWIN, RRANGE=RRANGE, SMTVAL=SMTVAL, OUTPUT=OUTPUT, $ QUIET=QUIET ; INITIALISE A FEW VARIABLES ;------------------------------------------------------------------------------- if n_elements(rrange) ne 2 then rrange = [7, 26] y_vals = findgen(rrange[1]-rrange[0]+1) + rrange[0] winsiz = size(detwin) dummyv = '' ; SUPPRESS ANY INFORMATIONAL MESSAGES IF DESIRED ;------------------------------------------------------------------------------- quieto = !quiet if keyword_set(quiet) then !quiet = 1 ; CREATE THE STRUCT CONTAINING THE FIT PARAMETERS ;------------------------------------------------------------------------------- tmpvar = fltarr(winsiz[1]) params = {factor:tmpvar, factorerr:tmpvar, positi:tmpvar, positierr:tmpvar, $ fwhmax:tmpvar, fwhmaxerr:tmpvar, chisqa:tmpvar, $ errflg:bytarr(winsiz[1])} ; OPTIONALLY SMOOTH THE DATA ;------------------------------------------------------------------------------- if n_elements(smtval) eq 2 then detwin = smooth(detwin, smtval, /edge_trunc) ; in future could use the following starting with IDL8.0: ;if n_elements(smtval) eq 2 then detwin = gauss_smooth(detwin, smtval, /edge_tr) ; CALCULATE SCALING PARAMETERS FOR PLOTTING AND SET PLOT ENVIRONMENT ;------------------------------------------------------------------------------- if keyword_set(output) then begin subset = detwin[*,rrange[0]:rrange[1]] sorted = sort(subset) yrange = subset[sorted[[0.02, 1.00] *n_elements(sorted)]] d_orig = !d.name p_orig = !p if !version.os_family eq 'Windows' then set_plot, 'WIN' $ else set_plot, 'X' window, /free, xsize=600, ysize=450 device, decompose=0 !p.charsize=1.25 endif ; GUESS THE ROW WHERE THE SPECTRUM IS LOCATED ;------------------------------------------------------------------------------- ; This does not work well if the spectrum is very curved, e.g. for GRISM!! tmpvar = max(total(detwin[*,rrange[0]:rrange[1]], 1), maxpos) ; Exteremely smooth the spectrum and then get position for each column smtwin = smooth(detwin, [50,0], /edge_trunc) ; LOOP OVER ALL COLUMNS I.E. WAVELENGTH BINS ;------------------------------------------------------------------------------- for i=0,winsiz[1]-1 do begin ; FIND THE MAXIMUM VALUE OF THE COLUMN (HOPING IT IS THE SPECTRUM) ;----------------------------------------------------------------------- ; This is the simplest guess of the amplitude of the Gaussian maxcol = max(detwin[i,rrange[0]:rrange[1]]) ; Use the extremely smoothed version to guess amplitude and position ; maxcol = max(smtwin[i,rrange[0]:rrange[1]], maxpos) ; FIT A GAUSSIAN TO THE COLUMN ;----------------------------------------------------------------------- tmpvar = gaussfit(y_vals, detwin[i,rrange[0]:rrange[1]], partmp, $ estimates = [maxcol, maxpos+rrange[0], 3, 0, 0], $ chisq=chisqa, nterms=5, sigma=errtmp) ; MAKE DIAGNOSTIC PLOTS FOR DEBUGGING OR CHECKING ;----------------------------------------------------------------------- if keyword_set(output) then begin plot, y_vals, detwin[i,rrange[0]:rrange[1]], yrange=yrange, $ xtitle='pixels on detector', ytitle='ADU', $ title='Column ' + string(i) oplot, y_vals, tmpvar, color=210 wait, 0.25 ;if i lt 50 then read, dummyv ;print, partmp endif ; COPY THE PARAMETERS INTO THE RESULT STRUCT ;----------------------------------------------------------------------- params.factor[i] = partmp[0] params.factorerr[i] = errtmp[0] params.positi[i] = partmp[1] params.positierr[i] = errtmp[1] params.fwhmax[i] = abs(partmp[2]) * 2 * sqrt(2*alog(2)) params.fwhmaxerr[i] = errtmp[2] * 2 * sqrt(2*alog(2)) params.chisqa[i] = chisqa ; TRY TO FIND BAD FITS AND FLAG THE DATA ACCORDINGLY ;----------------------------------------------------------------------- badfit = ((partmp[0] lt 0.0) or (params.fwhmax[i] lt 1.0) or $ (partmp[1] lt rrange[0]+2) or (partmp[1] gt rrange[1]-2)) if badfit then params.errflg[i] = 1B endfor ; RESTORE THE ORIGINAL PLOT SETTINGS ;------------------------------------------------------------------------------- if keyword_set(output) then begin set_plot, d_orig !p = p_orig endif ; RESET THE QUIET SYSTEM VARAIBLE TO ITS ORIGINAL VALUE ;------------------------------------------------------------------------------- !quiet = quieto ; RETURN THE RESULT STRUCT ;------------------------------------------------------------------------------- return, params END ;+ ; NAME: ; MIDI_MSKFIT ; ; PURPOSE: ; Create a MIDI mask from image detector files. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; mskdat = midi_mskfit(infile [, msknam] [, factor=factor] ; [, additi=additi] [, smtval=smtval] [, /output] ; [, psfile=psfile] [, skymsk=skymsk] ; [, skydis=skydis] [, skywid=skywid] ; [, linpar=linpar]) ; ; REQUIRED INPUTS: ; infile input FITS file with image detector data or IDL struct ; ; OPTIONAL INPUTS: ; msknam name of the mask file for saving the resulting mask ; ; KEYWORDS: ; factor a factor for enlarging the FWHM ; additi a constant for enlarging the FWHM ; smtval two element array with the smoothing values in x and y ; output plot the traces and FWHM of the mask ; psfile name of the PS file for PS output ; skydis distance of the sky band from the trace ; skywid width of the sky band in pixels ; linpar parameter for linear weighting in wavelength direction ; yshift shift the mask by the specified pixels in y-direction ; ; OUTPUTS: ; mskdat structure containing the object mask and its parameters ; skymsk structure containing the sky mask and its parameters ; ; DESCRIPTION AND EXAMPLE: ; This function takes either a name of an imaging detector FITS file ; or an IDL structure with imaging detector data as an input for the ; creation of a MIDI mask. If the data is complex (e.g. a fringeimage ; from oirAverageVisImage) then only the real part is considered, The ; input can have an arbitrary number of detector windows (e.g. for 2 ; for HIGH_SENSE or 4 for SCI_PHOT data). Multiple frames in the input ; data are interpreted as different beams (e.g. the two beams for A ; and B photometry). For each of the detector windows, the trace (i.e. ; the peak position as a function of column/wavelength) and the FWHM ; of the spectrum is determined by fitting a 1D Gaussian distribution ; to each column of the window. This is done for every beam seperately ; as well as for the average image of all beams. The trace and the FWHM ; are then parametrised as a function of column number by a polynomial. ; The final trace of the mask is defined as the average of the traces ; of the individual beams. The FWHM of the mask is defined as the ; maximum of the full width at half maxima of the different beams, ; multiplied by FACTOR, plus half of the distance between the traces ; and plus a constant ADDITI. The mask is then calculated as a modified ; Gaussian distribution ; z=exp( -0.5 * abs((y - trace) / FWHM*sqrt(2*alog(2)))^EXPONE ) ; in order to obtain a more bell shaped profile. The default value of ; EXPONE is 3. The result returned by this function is a structure of ; the same format as the input data extended by fields containing the ; parameters used for the creation of the mask. Additionally, a sky ; mask is calculated. The two sky bands are located symmetrically below ; and above the trace at a distance of SKYDIS * FWHM and the width of ; the sky bands is governed by the parameter SKYWID. If the keyword ; OUTPUT is set then a plot with the traces and the FWHM as well as ; comparisons between the final masks and the data are created. The ; masks can be written to FITS files by specifying a base for the file- ; names in MSKNAM. The parameters of the mask are stored in the FITS ; header as well as in the FITS table. A relatively important and ; currently hard-coded parameter is the range of rows used for fitting ; the spectrum. When fitting very faint spectra, then an optimal range ; of rows specified by this parameter is crucial to obtain a good fit, ; for bright sources (e.g. calibrators) this parameter is not so ; important as long as the spectrum falls within the range of rows. ; The function is simply called by: ; result = midi_mskfit(phtfil) ; Possible impürovements: (1) Add possibility for a providing a user ; defined weighting in wavelength direction (instead of LINPAR); ; (2) possibility to calculate a mask for given parameters, ; i.e. separating the mask fitting from the mask creating part. ; ; CALLED BY: ; none ; ; CALLING: ; mod_struct ; midi_wavaxi ; midi_gaudis ; ; MODIFICATION HISTORY: ; 2008-10-13 Written by Konrad R. Tristram ; 2008-11-19 K. Tristram: Added the possibility to smooth the windows. ; 2009-04-15 K. Tristram: Corrected error in handling of parameter RRANGE. ; 2009-04-26 K. Tristram: Added keyword ADDITI with default value of 1 ; pixel and changed the determination of the FWHM. ; 2009-12-21 K. Tristram: Added creation of a sky mask. ; 2010-07-16 K. Tristram: Added linear weighting in wavelength direction. ; 2010-11-05 K. Tristram: Major rewrite, now handles mutiple windows but ; only accepts reduced images. ; 2011-06-10 K. Tristram: Remove bug when writing out the FITS file. ; 2012-08-30 K. Tristram: Added keyword YSHIFT. ; FUNCTION MIDI_MSKFIT, INFILE, MSKNAM, FACTOR=FACTOR, ADDITI=ADDITI, $ SMTVAL=SMTVAL, OUTPUT=OUTPUT, PSFILE=PSFILE, $ SKYMSK=SKYMSK, SKYDIS=SKYDIS, SKYWID=SKYWID, $ YSHIFT=YSHIFT, LINPAR=LINPAR ; INITIALISE A FEW VARIABLES ;------------------------------------------------------------------------------- expone = 3 ; this gives a nice bell shape, similar to the standard EWS mask rrange = [ 8, 28] ; default range of rows (i.e. y-values) used for the fitting ;rrange = [ 6, 16] ; range of rows for UT data in 2009 ;rrange = [ 8, 20] ; range of rows for AT data in 2009 ;rrange = [10, 20] if n_elements(factor) ne 1 then factor = 1.25 ; This should be better 1.5!!! if n_elements(additi) ne 1 then additi = 2.00 if n_elements(smtval) ne 2 then smtval = [0., 0.] if n_elements(skydis) ne 1 then skydis = 1.0 if n_elements(skywid) ne 1 then skywid = 4.0 if n_elements(linpar) lt 1 then linpar = 0.0 ; DEFINE A VECTOR CONTAINING THE ATMOSPHERIC TRANSPARANCY ;------------------------------------------------------------------------------- transp = [0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, $ 0.000, 0.000, 0.000, 0.006, 0.015, 0.021, 0.030, 0.050, $ 0.080, 0.101, 0.105, 0.097, 0.088, 0.090, 0.099, 0.106, $ 0.114, 0.135, 0.180, 0.250, 0.353, 0.484, 0.601, 0.689, $ 0.754, 0.804, 0.842, 0.870, 0.888, 0.898, 0.907, 0.913, $ 0.914, 0.915, 0.922, 0.926, 0.928, 0.929, 0.917, 0.901, $ 0.895, 0.904, 0.925, 0.944, 0.955, 0.963, 0.969, 0.975, $ 0.980, 0.983, 0.982, 0.983, 0.987, 0.989, 0.989, 0.989, $ 0.988, 0.985, 0.979, 0.979, 0.987, 0.993, 0.995, 0.994, $ 0.994, 0.995, 0.995, 0.994, 0.993, 0.995, 0.999, 1.000, $ 0.999, 0.998, 0.999, 0.999, 0.998, 0.999, 1.000, 0.999, $ 0.996, 0.992, 0.991, 0.994, 0.996, 0.994, 0.988, 0.982, $ 0.975, 0.958, 0.924, 0.869, 0.793, 0.704, 0.614, 0.541, $ 0.488, 0.453, 0.447, 0.444, 0.429, 0.503, 0.685, 0.851, $ 0.930, 0.954, 0.962, 0.963, 0.960, 0.955, 0.949, 0.944, $ 0.944, 0.944, 0.941, 0.934, 0.934, 0.942, 0.943, 0.925, $ 0.901, 0.861, 0.771, 0.627, 0.469, 0.353, 0.285, 0.251, $ 0.247, 0.249, 0.226, 0.170, 0.103, 0.055, 0.024, 0.007, $ 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, $ 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.009, 0.025, $ 0.033, 0.023, 0.005, 0.000, 0.000, 0.000, 0.000, 0.000, $ 0.000, 0.000, 0.000] ; READ THE DATA OR ASSUME IT IS ALREADY A STRUCT OF THE CORRECT TYPE ;------------------------------------------------------------------------------- if size(infile, /type) eq 7 then begin indata = oirgetdata(infile) refile = (strsplit(infile[0], /extract))[0] endif else indata = infile numele = n_elements(indata) ; FIND HOW MANY AND WHERE THE DATA WINDOWS ARE ;------------------------------------------------------------------------------- tagnam = tag_names(indata) tagnum = n_tags(indata) winidx = where(strmatch(tagnam, 'DATA*')) if winidx[0] lt 0 then begin message, 'No detector data fields detected in the data. Exiting!', /cont return, -1 endif winidx = winidx[sort(tagnam[winidx])] winnum = n_elements(winidx) ; IF DATA IS PSEUDOCOMPLEX (E.G. FROM FRINGEIMAGE) THEN MAKE IT FLOAT ;------------------------------------------------------------------------------- winsiz = size(indata[0].data1) if (winsiz[1] eq 342) or (winsiz[1] eq 522) then begin ; CREATE A NEW STRUCTURE AND MODIFY THE SIZE OF THE DATA FIELDS ;----------------------------------------------------------------------- tmpvar = indata[0] for i=0,winnum-1 do tmpvar = mod_struct(tmpvar, tagnam[winidx[i]], $ fltarr(winsiz[1]/2, winsiz[2])) tmpvar = replicate(tmpvar, numele) ; CALCULATE THE REAL PART OF THE DATA AND COPY IT INTO THE NEW STRUCTURE ;----------------------------------------------------------------------- for i=0,winnum-1 do tmpvar.(winidx[i]) = (-1)^i * $ transpose(float(pseudocomplex(indata.(winidx[i])))) ; MOVE THE RESULT TO THE INPUT DATA STRUCTURE ;----------------------------------------------------------------------- indata = temporary(tmpvar) endif ; FIND OUT THE DATA TYPE (GRISM OR PRISM) ;------------------------------------------------------------------------------- winsiz = size(indata[0].data1) case winsiz[1] of 151: begin disper = 'prismold' colmin = 05 colmax = 105 colexa = [ 20, 90] yrange = [ 12, 32] end 171: begin disper = 'prism' colmin = 30 colmax = 140 colexa = [ 50, 120] yrange = [ 8, 22] end 261: begin disper = 'grism' colmin = 20 colmax = 250 colexa = [ 30, 200] yrange = [ 5, 25] end else: begin print, 'Unknown dispersion type or window type! Exiting.' return, -1 end endcase ; DETERMINE THE REFERENCE FILE IF NOT ALREADY DONE ;------------------------------------------------------------------------------- if winnum ge 3 then tmpvar='smask' else tmpvar='hmask' if n_elements(refile) lt 1 then refile = getenv(disper+tmpvar) ; GET THE WAVELENGTH SCALING ;------------------------------------------------------------------------------- wavele = oirgetwavelength(refile) ; CALCULATE A FEW ARRAYS DEPENDING ON THE DATA TYPE ;------------------------------------------------------------------------------- coln_1 = findgen(colmax-colmin+1)+colmin coln_2 = findgen(winsiz[1]) rownum = findgen(winsiz[2]) ; CREATE THE STRUCT HOLDING THE RESULTING MASKS ;------------------------------------------------------------------------------- newmsk = indata[0] for i=1,winnum do newmsk = create_struct(newmsk, $ 'TRACE'+strtrim(i, 2), fltarr(3), 'FWHM_'+strtrim(i, 2), fltarr(3)) ; CHECK IF INPUT IS INTEGER, THEN CHANGE THE MASK TO FLOAT ;------------------------------------------------------------------------------- for i=0,winnum-1 do if size(newmsk.(winidx[i]), /type) eq 2 then $ newmsk = mod_struct(newmsk, tagnam[winidx[i]], 1.0*newmsk.(winidx[i])) ; NOW ALSO CREATE THE SKY MASK, WHICH IS IDENTICAL TO THE MASK IN STRUCTURE ;------------------------------------------------------------------------------- skymsk = newmsk ; CALCULATE WEIGHTS DEPENDING LINEARLY ON WAVELENGTH ;------------------------------------------------------------------------------- if linpar[0] ge 0. then slopef = 1. - linpar[0]*(max(wavele)-wavele) < 1. > 0. $ else slopef = 1. - linpar[0]*(min(wavele)-wavele) < 1. > 0. ; SET PLOTTING ENVIRONMENT DEPENDING ON THE DEVICE ;------------------------------------------------------------------------------- if keyword_set(output) then begin d_orig = !d.name p_orig = !p if keyword_set(psfile) then begin set_plot, 'PS' device, file=psfile+'.ps', bits_per_pixel=8, /landscape, $ xoffset=1.0, yoffset=28.3, xsize=27., ysize=19., /color !p.charsize=1.5 colors = [60,100] endif else begin if !version.os_family eq 'Windows' then set_plot, 'WIN' $ else set_plot, 'X' window, /free, xsize=1200, ysize=800 device, decompose=0 !p.charsize=2.5 colors = [100, 60] endelse !p.multi = [0, 3, 2] loadct, 39, /silent endif ; LOOP OVER THE DETECTOR WINDOWS ;------------------------------------------------------------------------------- for i=0,winnum-1 do begin ; GET THE TAG INDEX FOR THE DATA OF THE CURRENT WINDOW ;----------------------------------------------------------------------- curidx = winidx[i] ; FIT GAUSSIANS COLUMNWISE TO THE COMBINATION OF THE BEAMS ;----------------------------------------------------------------------- if numele gt 1 then tmpvar = total(indata.(curidx), 3) $ else tmpvar = indata.(curidx) params = midi_gaudis(tmpvar, smtval=smtval, rrange=rrange, /quiet) ; NOW FIT GAUSSIANS COLUMNWISE TO ALL BEAMS INDIVIDUALLY ;----------------------------------------------------------------------- params = replicate(params, numele+1) for j=0,numele-1 do params[j] = midi_gaudis(indata[j].(curidx), $ smtval=smtval, rrange=rrange, /quiet) ; SET THE ERRORS FOR BAD FITS TO INFINITY ;----------------------------------------------------------------------- traces = replicate({p:fltarr(3), v:fltarr(winsiz[1])}, numele+1) fwhmax = replicate({p:fltarr(3), v:fltarr(winsiz[1])}, numele+2) for j=0,numele do begin ; SET THE ERRORS FOR BAD FITS TO INFINITY ;--------------------------------------------------------------- badfit = where(params[j].errflg eq 1B) > 0 params[j].positierr[badfit] = !values.f_infinity params[j].fwhmaxerr[badfit] = !values.f_infinity ; FIT THE TRACES BY POLYNOMIALS ;--------------------------------------------------------------- traces[j].p = poly_fit(coln_1, params[j].positi[colmin:colmax],$ 2, measure=params[j].positierr[colmin:colmax]) traces[j].v = poly(coln_2, traces[j].p) ; FIT THE FWHMS BY POLYNOMIALS ;--------------------------------------------------------------- fwhmax[j].p = poly_fit(coln_1, params[j].fwhmax[colmin:colmax],$ 2, measure=params[j].fwhmaxerr[colmin:colmax]) fwhmax[j].v = poly(coln_2, fwhmax[j].p) endfor ; FIT THE SEPARATION OF THE TRACES BY A POLYNOMIAL ;----------------------------------------------------------------------- ; The following fits the separation of the traces: tmpvar = 0.5* (max(traces.v, dim=2) - min(traces.v, dim=2)) ; The following directly fits the separation of the data: ; tmpvar = 0.5* (max(params.positi, maxidx, dim=2) - $ ; min(params.positi, minidx, dim=2)) separa = poly_fit(coln_1, tmpvar[colmin:colmax], 2) ; COPY THE AVERAGED TRACE PARAMETERS INTO THE RESULTING MASK STRUCT ;----------------------------------------------------------------------- newmsk.(tagnum+2*i+0) = total(reform(traces[0:numele-1].p, $ 3, numele), 2) / numele ; APPLY AN ADDITIONAL SHIFT TO THE MASK IN Y-DIRECTION ;----------------------------------------------------------------------- if keyword_set(yshift) then newmsk.(tagnum+2*i+0)[0] += yshift ; PRINT OUT THE COEFFICIENTS AND MAKE A DIAGNOSTIC PLOT FOR THE TRACE ;----------------------------------------------------------------------- if keyword_set(output) then begin print, format='("Window", I2, ": Trace y =", F8.4, ' + $ 'F+11.6, " x", E+13.4, " x²")', i+1, $ newmsk.(tagnum+2*i+0) !p.multi = [3*(winnum-i), 3, winnum] plot, [0,0], xstyle=9, xrange=[0,winsiz[1]], xmargin=[6, 2], $ ystyle=1, yrange=yrange, /nodata, $ xtitle='x [pixels]', ytitle='y [pixels]', $ title='Trace for window ' + strtrim(i+1,1) midi_wavaxi, disper, xtickl=0.025 oplot, [colexa[0], colexa[0]], [-100,100], linestyle=1 oplot, [colexa[1], colexa[1]], [-100,100], linestyle=1 for j=0,numele-1 do oplot, params[j].positi, color=colors[0] oplot, params[numele].positi, color=colors[1] for j=0,numele-1 do oplot, traces[j].v, color=200 oplot, poly(coln_2, newmsk.(tagnum+2*i+0)), thick=2, color=250 oplot, poly(coln_2, separa)+yrange[0], linestyle=1, color=250 xyouts, winsiz(1)/2., 0.083*yrange[0]+ 0.916*yrange[1], $ string(format='(" y!D", I1, "!N =", F6.2, ' + $ 'F+7.3, " x", E+10.2, " x!U2!N")', i+1, $ newmsk.(tagnum+2*i+0)), charsize=0.4*!p.charsize, $ align=0.5 endif ; FIT THE ENVELOPE OF THE FWHM BY A POLYNOMIAL ;----------------------------------------------------------------------- tmpvar = max(params.fwhmax, maxidx, dim=2) fwhmax[numele+1].p = poly_fit(coln_1, tmpvar[colmin:colmax], 1, $ measure=((params.fwhmaxerr)[maxidx])[colmin:colmax]) fwhmax[numele+1].v = poly(coln_2, fwhmax[numele+1].p) ; COPY THE TRACE INTO THE RESULTING MASK STRUCT ;----------------------------------------------------------------------- newmsk.(tagnum+2*i+1) = factor * fwhmax[numele+1].p + 2.0 * separa + $ [additi, 0.0, 0.0] ; PRINT OUT THE COEFFICIENTS AND MAKE A DIAGNOSTIC PLOT FOR THE FWHM ;----------------------------------------------------------------------- if keyword_set(output) then begin print, format='(" FWHM =", F8.4, F+11.6, ' + $ '" x", E+13.4, " x²")', newmsk.(tagnum+2*i+1) plot, [0,0], xstyle=9, xrange=[0,winsiz[1]], xmargin=[6, 2], $ ystyle=1, yrange=[0,14], /nodata, $ xtitle='x [pixels]', ytitle='FWHM y [pixels]', $ title='FWHM for window ' + strtrim(i+1,1) midi_wavaxi, disper, xtickl=0.025 oplot, [colexa[0], colexa[0]], [-100,100], linestyle=1 oplot, [colexa[1], colexa[1]], [-100,100], linestyle=1 for j=0,numele-1 do oplot, params[j].fwhmax, color=colors[0] oplot, tmpvar, color=colors[1] for j=0,numele-1 do oplot, fwhmax[j].v, color=200 oplot, fwhmax[numele+1].v, color=250 oplot, poly(coln_2, newmsk.(tagnum+2*i+1)), thick=2, color=250 xyouts, winsiz(1)/2., 0.083*0.0+ 0.916*14.0, $ string(format='(" FWHM!D", I1, "!N =", F6.2, ' + $ 'F+7.3, " x", E+10.2, " x!U2!N")', i+1, $ newmsk.(tagnum+2*i+1)), charsize=0.4*!p.charsize, $ align=0.5 endif ; COPY THE MASK PARAMETERS FOR THIS WINDOW ALSO INTO THE SKY MASK ;----------------------------------------------------------------------- for j=0,1 do skymsk.(tagnum+2*i+j) = newmsk.(tagnum+2*i+j) ; CALCULATE THE MASK FOR THE SPECTRUM AND THE SKY ;----------------------------------------------------------------------- ; an easy way to limit the mask in wavelength direction is to only ; calculate it for a certain range of columns ; for j=25,140 do begin for j=0,winsiz[1]-1 do begin ; FIRST CALCULATE THE MASK FOR THE SPECTRUM ;--------------------------------------------------------------- tmpvar = (abs((rownum - poly(j, newmsk.(tagnum+2*i+0))) / $ poly(j, newmsk.(tagnum+2*i+1)) * $ 2*sqrt(2*alog(2))))^expone newmsk.(curidx)[j,*] = exp(-0.5*tmpvar) * slopef[j]; * transp[j] ; EITHER USE VARIABLE OR FIXED SKY DISTANCE ;--------------------------------------------------------------- if 1 then skysep = skydis * poly(j, newmsk.(tagnum+2*i+1)) $ else skysep = skydis ; NOW CALCULATE THE MASK FOR THE SKY BANDS ;--------------------------------------------------------------- tmpvar = exp(-0.5 * (2.0 / skywid * $ abs((rownum-poly(j, newmsk.(tagnum+2*i)) - skysep)))^10.) tmpvar += exp(-0.5 * (2.0 / skywid * $ abs((rownum-poly(j, newmsk.(tagnum+2*i)) + skysep)))^10.) skymsk.(curidx)[j,*] = tmpvar endfor newmsk.(curidx) = (1.1 * newmsk.(curidx) - 0.05) < 1.0 > 0.0 skymsk.(curidx) = (2.0 * skymsk.(curidx) - 1.00) < 1.0 > 0.0 ; PLOT COMPARISONS BETWEEN THE PROFILE OF THE SPECTRUM AND THE MASK ;----------------------------------------------------------------------- if keyword_set(output) then begin !p.multi = [(winnum-i-1)*6+4, 3, 2*winnum] plot, [0,35], [-0.2,1.19], xstyle=1, xmargin=[6,2], $ ystyle=1, ymargin=[1,2], xtickname=replicate(' ', 30), $ xticklen=0.04, title='Cuts at columns ' + $ strtrim(colexa[0],1) + ' and ' + strtrim(colexa[1],1), $ /nodata for j=0B,numele-1 do begin tmpvar = float(indata[j].(curidx)[colexa[0],*]) oplot, tmpvar / max(tmpvar[rrange[0]:rrange[1]], $ maxidx), color=colors[0] if numele ge 2 then xyouts, maxidx+rrange[0]-0.5, $ 1.05, string(j+65B), color=colors[0], $ charsize=0.75*!p.charsize endfor oplot, newmsk.(curidx)[ colexa[0],*], color=250 oplot, skymsk.(curidx)[ colexa[0],*], color=200 !p.multi = [(winnum-i-1)*6+1, 3, 2*winnum] plot, [0,35], [-0.2,1.19], xstyle=1, xmargin=[6,2], $ ystyle=1, ymargin=[4,-1], xticklen=0.04, $ xtitle= 'y [pixels]', /nodata for j=0B,numele-1 do begin tmpvar = float(indata[j].(curidx)[colexa[1],*]) oplot, tmpvar / max(tmpvar[rrange[0]:rrange[1]], $ maxidx), color=colors[0] if numele ge 2 then xyouts, maxidx+rrange[0]-0.5, $ 1.05, string(j+65B), color=colors[0], $ charsize=0.75*!p.charsize endfor oplot, newmsk.(curidx)[ colexa[1],*], color=250 oplot, skymsk.(curidx)[ colexa[1],*], color=200 endif endfor ; CLOSE THE DEVICE IF CREATING A PS FILE ;------------------------------------------------------------------------------- if keyword_set(psfile) then device, /close ; RESTORE THE ORIGINAL PLOT SETTINGS ;------------------------------------------------------------------------------- if keyword_set(output) then begin !p = p_orig set_plot, d_orig endif ; SAVE THE MASK AS A FITS FILE IF A NAME HAS BEEN GIVEN ;------------------------------------------------------------------------------- if size(msknam, /type) eq 7 then begin ; OLD AND SIMPLE OUTPUT OF THE MASK INTO A FITS FILE ;----------------------------------------------------------------------- ; oirMakeMask, refile, msknam+'-spc1.fits', newmsk ; oirMakeMask, refile, msknam+'-sky1.fits', skymsk ; SET UP THE MAIN HEADER OF THE MASKFILES, ADJUST NAXIS ;----------------------------------------------------------------------- refdet = obj_new('imagedetector', refile) header = (refdet->priHead())->clone() refdet->close header->addpar, 'OBJECT', 'Object mask', comment=' Original target' header->addpar, 'ORIGIN', 'MIDI_MSKFIT', $ commen=' MIDI_MSKFIT mask fitting routine v2010-11-08' reftab = refdet->table() for i=0,winnum-1 do (*reftab)[i].naxis[0] = winsiz[1] ; ADD THE PARAMETERS OF THE FUNCTION CALL TO THE MAIN HEADER ;----------------------------------------------------------------------- comstr = ' MIDI_MSKFIT: ' header->addpar, 'M_FACTOR', factor, com=comstr+'broadening factor' header->addpar, 'M_ADDITI', additi, com=comstr+'additional width' header->addpar, 'M_SMOOTX', smtval[0], com=comstr+'smoothing width in x' header->addpar, 'M_SMOOTY', smtval[1], com=comstr+'smoothing width in y' header->addpar, 'M_SKYDIS', skydis, com=comstr+'sky distance' header->addpar, 'M_SKYWID', skywid, com=comstr+'sky band width' keywrd = ['(constant)', '(linear)', '(quadratic)'] ; ADD THE PARAMETERS OF THE MASK FITS TO THE MAIN HEADER ;----------------------------------------------------------------------- for i=0,winnum-1 do begin tmpstr = ' window '+strtrim(i+1,2)+' ' for j=0B,2B do header->addpar, $ 'M_T'+strtrim(i+1,2)+'COE'+strtrim(j+1,2), $ (newmsk.(tagnum+2*i+0))[j], comment=comstr+'Trace' + $ tmpstr+'coeff '+strtrim(j+1,2)+' '+keywrd[j] for j=0B,2B do header->addpar, $ 'M_W'+strtrim(i+1,2)+'COE'+strtrim(j+1,2), $ (newmsk.(tagnum+2*i+1))[j], comment=comstr+'Width' + $ tmpstr+'coeff '+strtrim(j+1,2)+' '+keywrd[j] endfor ; WRITE OUT THE DATA INTO FITS FILES FOR UNCHANGED FITS TABLE STRUCTURE ;----------------------------------------------------------------------- ; This does not work if the input is complex, i.e. a fringe image! ; header->addpar, 'OBJECT', 'Object mask', comment=' Original target' ; oirNewData, refile, msknam+'-spc2.fits', newmsk, newHead=header ; header->addpar, 'OBJECT', 'Sky mask', com=' Original target' ; oirNewData, refile, msknam+'-sky2.fits', skymsk, newHead=header ; ADD REQUIRED AND ADDITIONAL KEYWORDS TO THE HEADER OF THE FITS TABLE ;----------------------------------------------------------------------- output = obj_new('fitstable', newmsk, extname='IMAGING_DATA') output->addpar, 'ORIGIN', 'MIDI_MSKFIT', $ comment=' MIDI_MSKFIT mask fitting routine v2010-11-08' output->addpar, 'INSTRUME', 'MIDI ', comment=' Instrument used' output->addpar, 'NREGION', winnum, comment=' The number of regions ' + $ 'used in the current mode' output->addpar, 'MAXTEL', 2, comment=' The maximum number of ' + $ 'telescopes contributing p' output->addpar, 'MAXINS', 1, comment=' Number of entries in ' + $ 'INS_TRAIN table' output->addpar,'REVISION', 1, comment=' Revision number of the ' + $ 'table definitio' ; WRITE OUT THE MASK FOR THE SPECTRUM ;----------------------------------------------------------------------- header->addpar, 'OBJECT', 'Object mask', comment=' Original target' refdet->newfile, msknam+'.srcmsk.fits', priHead=header->clone() output->appendtofile, refdet output->writerows, newmsk output->close ; WRITE OUT THE MASK FOR THE SKY BANDS ;----------------------------------------------------------------------- header->addpar, 'OBJECT', 'Sky mask', comment=' Original target' refdet->newfile, msknam+'.skymsk.fits', priHead=header output->appendtofile, refdet output->writerows, skymsk output->close ; CLEAN UP THE OBJECTS ;----------------------------------------------------------------------- refdet->close obj_destroy, refdet obj_destroy, output obj_destroy, header endif ; RETURN THE RESULT MASK ;------------------------------------------------------------------------------- return, newmsk END ;+ ; NAME: ; MIDI_FLAGOK ; ; PURPOSE: ; Determine which frames where flagged as good in a fringe track. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; flagok = midi_flagok(indata [, frames]) ; ; REQUIRED INPUTS: ; indata EWS tag or structure containing flagging information ; ; OPTIONAL INPUTS: ; frames structure containg the frame data ; ; KEYWORDS: ; none ; ; OUTPUTS: ; result ; ; DESCRIPTION AND EXAMPLE: ; This function determines which frames of a fringe track where flagged ; as good by oirAutoFlag and returns an array with all good frames ; assigned the value 1, all others the value 0. The input data INDATA ; to the routine can either be an IDL string or a structure. If it is ; an IDL string, INDATA is interpreted as an EWS tag and the flagging ; as well as the frame data are directly loaded from the respective ; files (if present): tag.flag.fits and tag.fringes.fits. If INDATA is ; an IDL structure, it is interpreted as the structure containg the ; flagging data. In this case also the frame data has to be directly ; provided to the function using the second input parameter FRAMES. ; FRAMES can then also be the group delay data from oirGroupDelay with ; a reduced number of elements when using the slope fitting routines. ; The function is called by: ; midi_flagok, ewstag ; ; CALLED BY: ; midi_pltphi ; midi_pltopd ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2011-06-15 Created from parts of MIDI_PLTOPD by Konrad R. Tristram ; 2011-08-03 K. Tristram: Speeded up programme by copying variables ; 2011-11-23 K. Tristram: Add a check if FITSTABLE-object is valid ; 2012-04-03 K. Tristram: Changed counters in loops variables to LONG. ; FUNCTION MIDI_FLAGOK, INDATA, FRAMES ; CHECK WHAT KIND OF INPUT DATA WE ARE DEALING WITH ;------------------------------------------------------------------------------- case size(indata, /type) of ; IF IT IS A STRING THEN IT IS AN EWS TAG THEN READ THE DATA ;----------------------------------------------------------------------- 7: if file_test(indata+'.flag.fits') and $ file_test(indata+'.fringes.fits') then begin tmpvar = obj_new('fitstable', indata+'.flag.fits', $ extname='FLAG') if obj_valid(tmpvar) then begin flgdat = tmpvar->readrows() obj_destroy, tmpvar endif else begin obj_destroy, tmpvar message, 'Data seems corrupted. Exiting!', /cont return, -1 endelse ftimes = (oirgetdata(indata+'.fringes.fits')).time endif else begin message, 'Necessary files not found. Exiting!', /cont return, -1 endelse ; IF IT IS A STRUCT THEN COPY IT INTO THE RIGHT VARIABLES ;----------------------------------------------------------------------- 8: if size(frames, /type) eq 8 then begin flgdat = indata ftimes = frames.time endif else begin message, 'Also have to specify GDELAY. Exiting!', /cont return, -1 endelse else: begin message, 'INPUT has incorrect data format. Exiting!', /cont return, -1 end endcase ; DETERMINE THE NUMBER OF GROUP DELAY ESTIMATES ;------------------------------------------------------------------------------- number = n_elements(ftimes) ; DETERMINE ALL ELEMENTS WHICH ARE FLAGGED AS GOOD DATA ;------------------------------------------------------------------------------- flagok = bytarr(number) + 1B flgdur = flgdat.timerang[1] - flgdat.timerang[0] for i=0L,n_elements(flgdat)-1 do if flgdur[i] gt 0.0 then begin flgsub = where((ftimes ge flgdat[i].timerang[0]) and $ (ftimes le flgdat[i].timerang[1])) if flgsub[0] gt -1 then flagok[flgsub < (number-1)] = 0 endif ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, flagok END ;+ ; NAME: ; MIDI_UNWRAP ; ; PURPOSE: ; Unwrap phases assuming a smooth phase change ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; newphi = midi_unwrap(oldphi [, smooth] [, maxval=40]) ; ; REQUIRED INPUTS: ; oldphi array containing the input phases (in degrees) ; ; OPTIONAL INPUTS: ; smooth smoothing width to use for finding the phase values ; ; KEYWORDS: ; maxval maximum allowed jump before assuming a phase wrap ; ; OUTPUTS: ; phases ; ; DESCRIPTION AND EXAMPLE: ; This function unwraps phases in an "intelligent" way, i.e. not ; getting confused by outliers. The value of every phase in the ; array is compared to the previous phases. If a difference larger ; than MAXVAL is detected, a phase wrap is assumed and the following ; phases are corrected for this wrap. In order to not follow ; individual outliers, the comparison is carried out to the mean of ; a number of previous phases. The number of phases used in this ; mean is specified by the keyword SMOOTH. Phases are assumed to be ; in degrees. ; The function is called by: ; newphi = midi_unwrap(oldphi) ; ; CALLED BY: ; midi_getphi ; midi_pltphi ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2011-12-16 Written by Konrad R. Tristram ; 2012-04-03 K. Tristram: Changed counters in loops variables to LONG. ; FUNCTION MIDI_UNWRAP, OLDPHI, SMOOTH, MAXVAL=MAXVAL ; DEFINE THE DEFAULT VALUES FOR THE SMOOTHING AND MAXIMUM ALLOWED JUMPS ;------------------------------------------------------------------------------- if n_elements(smooth) ne 1 then smooth = 40 $ ; one piezo scan with MIDI else smooth = smooth > 1 if n_elements(maxval) ne 1 then maxval = 200 ; a little more than 180° ; COPY THE INPUT ARRAY TO THE OUTPUT ARRAY ;------------------------------------------------------------------------------- newphi = oldphi ; LOOP OVER EVERY ELEMENT OF THE ARRAY ;------------------------------------------------------------------------------- for i=1L,n_elements(newphi)-1 do begin ; CALCULATE THE DIFFERENCE TO THE PREVIOUS PHASES ;----------------------------------------------------------------------- differ = newphi[i] - mean(newphi[(i-smooth)>0:i-1]) ; CHECK IF A WRAP OCURRES AND CORRECT FOR IT ;----------------------------------------------------------------------- if differ gt +maxval then newphi[i:*] -= 360. if differ lt -maxval then newphi[i:*] += 360. endfor ; SHIFT THE ENTIRE PHASE ARRAY CLOSE TO PHASE 0 ;------------------------------------------------------------------------------- sorted = sort(newphi) midpos = 0.5 * total(newphi[sorted[[0.01, 0.99]*n_elements(newphi)]]) newphi -= ((fix(midpos) / 270) * 360) ; RETURN THE UNWRAPPED PHASES ;------------------------------------------------------------------------------- return, newphi END ;+ ; NAME: ; MIDI_GETPHI ; ; PURPOSE: ; Get the water vapour phases from a FITS file. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; phases = midi_getphi(filnam [, /unwrap]) ; ; REQUIRED INPUTS: ; filnam EWS tag or file name for the file containing the phases ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; unwrap unwrap the phases using MIDI_UNWRAP ; ; OUTPUTS: ; phases ; ; DESCRIPTION AND EXAMPLE: ; This function reads the water vapour phases from the FITS file ; produced as an output of oirRotateGroupDelay. The function can ; be either called with the full name of the FITS file or by simply ; providing the EWS tag. By setting the keyword UNWRAP, the phases ; will be unwrapped, i.e. jumps of 360° will be removed. ; The function is called by: ; phases = midi_getphi('/tmp/tmp.ungroupdelay.fits') ; ; CALLED BY: ; midi_pltphi ; ; CALLING: ; midi_unwrap ; ; MODIFICATION HISTORY: ; 2011-12-16 Written by Konrad R. Tristram ; 2012-07-19 K. Tristram: Catch erroneous time stamps ; FUNCTION MIDI_GETPHI, FILNAM, UNWRAP=UNWRAP ; DETERMINE IF ONLY EWSTAG WAS GIVEN OR PERHAPS A FULL FILENAME ;------------------------------------------------------------------------------- if strmatch(strmid(filnam, 4, 5, /reverse), '.fits') then begin tmpvar = strsplit(filnam, '.', /extract) ewstag = tmpvar[0] phifil = '.'+strjoin(tmpvar[1:*], '.') endif else begin ewstag = filnam phifil = '.ungroupdelay.fits' ; standard file for the phase data endelse ; CHECK IF THE INPUT FILE IS PRESENT ;------------------------------------------------------------------------------- if ~ file_test(ewstag+phifil) then begin message, 'File containing the phases does not exist. Exiting!', /cont return, -1 endif ; READ IN THE FITS TABLE CONTAINING THE DISPERSION PHASES ;------------------------------------------------------------------------------- tblobj = obj_new('fitstable', ewstag+phifil, extname='DISPERSION') if tblobj eq obj_new() then begin message, 'File does not contain a DISPERSION table. Exiting!', /cont obj_destroy, tblobj return, -1 endif phases = tblobj->readrows() obj_destroy, tblobj ; CATCH ERRONEOUS TIME STAMPS IN THE DISPERSION TABLE ;------------------------------------------------------------------------------- if (min(phases.time) lt 1.0E1) or (max(phases.time) gt 1.0E5) then begin phases.time = (oirgetdata(ewstag+phifil)).time endif ; CORRECT PHASE JUMPS WHEN THE PHASE WRAPS ;------------------------------------------------------------------------------- if keyword_set(unwrap) then phases.dispersion = midi_unwrap(phases.dispersion) ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, phases END ;+ ; NAME: ; MIDI_PLTPHI ; ; PURPOSE: ; Plot the water vapour phases for a MIDI fringe track. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; midi_pltphi, filnam [, xrange=xrange] [, yrange=yrange] ; [, noflag=noflag] [, titles=titles] ; [, symbol=symbol] [, refphi=refphi] ; [, smtval=smtval] [, pltpos=pltpos] ; [, noxlab=noxlab] [, psfile=psfile] ; ; REQUIRED INPUTS: ; filnam EWS tag or file name of the MIDI phase data ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; xrange range of the x axis (time in seconds) ; yrange range of the y axis (phase in degrees) ; noflag do not mark the flagged frames by a different colour ; titles string containing a custom title for the plot ; symbol mark the individual measurements with small points ; refphi reference phases to be plotted in red ; smtval smoothing value for the unwrapping of the phases ; pltpos manually specified device coordinates for the plot position ; noxlab suppress labeling of lower x axis ; psfile name of the PS file for PS output ; ; OUTPUTS: ; a plot of the water vapour phases ; ; DESCRIPTION AND EXAMPLE: ; This procedure loads the water vapour phases (from the dispersion ; table in *.ungroupdelay.fits) and the flagging information (from ; *.flag.fits) for an EWS tag specified in FILNAM in order to produce ; a plot of the phases as a function of time. The phases of frames ; flagged as "good" are overplotted in a different colour over the ; continuous phases, unless the keyword NOFLAG is set. The phases ; are plotted repeatedly with offsets of 360° in the background. By ; default, the plot is scaled to contain the entire fringe track and ; to range from -1000 to 1000 degrees. Custom ranges for the x and ; y axis can be specified with the XRANGE and YRANGE keywords. ; Reference phases can be provided using the keyword REFDEL, which ; must be an array of structures as obtained when using MIDI_GETPHI. ; In its simplest form, the procedure is called by: ; midi_pltphi, 'tag' ; ; CALLED BY: ; midi_pltres ; ; CALLING: ; midi_flagok ; midi_getphi ; midi_unwrap ; setplot ; ; MODIFICATION HISTORY: ; 2011-12-16 Written by Konrad R. Tristram ; 2012-06-04 K. Tristram: Added keywords PLTPOS and NOXLAB ; 2012-08-03 K. Tristram: Changed automatic calculation of XRANGE. ; PRO MIDI_PLTPHI, FILNAM, XRANGE=XRANGE, YRANGE=YRANGE, NOFLAG=NOFLAG, $ TITLES=TITLES, SYMBOL=SYMBOL, REFPHI=REFPHI, $ SMTVAL=SMTVAL, PLTPOS=PLTPOS, NOXLAB=NOXLAB, PSFILE=PSFILE ; DETERMINE IF ONLY EWSTAG WAS GIVEN OR PERHAPS A FULL FILENAME ;------------------------------------------------------------------------------- if strmatch(strmid(filnam, 4, 5, /reverse), '.fits') then begin tmpvar = strsplit(filnam, '.', /extract) ewstag = tmpvar[0] phifil = '.'+strjoin(tmpvar[1:*], '.') endif else begin ewstag = filnam phifil = '.ungroupdelay.fits' ; standard file for the phase data endelse ; GET THE PHASES FROM THE FILE UND UNWRAP THEM ;------------------------------------------------------------------------------- phases = midi_getphi(ewstag+phifil) if size(phases, /type) ne 8 then return unwrap = midi_unwrap(phases.dispersion, smtval) ; LOAD THE FLAGGING INFORMATION ;------------------------------------------------------------------------------- if (~ keyword_set(noflag)) then if file_test(ewstag+'.flag.fits') then begin tmpvar = obj_new('fitstable', ewstag + '.flag.fits', extname='FLAG') flgdat = tmpvar->readrows() obj_destroy, tmpvar message, 'Loaded water vapour phases from '+ewstag+phifil+'.', /continue message, 'Loaded flagging data from '+ewstag+'.flag.fits.', /conti endif else begin message, 'No flagging data found.', /continue noflag = 1B endelse ; IF POINTS ARE DESIRED THEN CREATE A POINT AS A USER SYMBOL ;------------------------------------------------------------------------------- if keyword_set(symbol) then begin usersym, cos(findgen(17) * (!pi*2/16.)), $ sin(findgen(17) * (!pi*2/16.)), /FILL pltsym = -8 endif else pltsym = 0 ; MAKE A TITLE FOR THE PLOT IF NONE HAS BEEN PROVIDED ;------------------------------------------------------------------------------- if ~ keyword_set(titles) then begin titles = strmid(oirGetKey(ewstag+phifil, 'DATE-OBS'), 0, 19) + $ ' - ' + oirGetKey(ewstag+phifil, 'DET NRTS MODE') + ': ' + $ oirGetKey(ewstag+phifil, 'OBS TARG NAME') endif ; CALCULATE THE TIMES AND SET DEFAULT VALUES FOR THE XRANGE AND YRANGE ;------------------------------------------------------------------------------- i_time = reform(phases.time-phases[0].time[0])*24*3600 if (n_elements(xrange) ne 2) then xrange=i_time[[0, n_elements(i_time)-1]] if (n_elements(yrange) ne 2) then yrange=[-1000,+1000] ; SET THE PLOTTING ENVIRONMENT IF NOT PROVIDED AND LOAD THE COLOUR TABLE ;------------------------------------------------------------------------------- if n_elements(pltpos) ne 4 then begin setplot, pltset, pltpos, /longer, psfile=psfile pltpos *= [1.0, 1.0, 1.0, 0.94] endif loadct, 0, /silent ; SET THE SYMBOL SIZE DEPENDING ON IF IT IS A PS PLOT OR NOT ;------------------------------------------------------------------------------- if !d.name eq 'PS' then symsiz=0.30 else symsiz=0.60 ; CREATE THE THE COORDINATE SYSTEM OF THE PLOT ;------------------------------------------------------------------------------- plot, xrange, yrange, yrange=yrange, xstyle=5, ystyle=5, position=pltpos, $ /device, /nodata, /noerase ; PLOT MULTIPLE PHASE TRACKS IN THE BACKGROUND ;------------------------------------------------------------------------------- if !d.name eq 'PS' then colour = [170, -1] else colour = [50, +1] for i=-2,2 do oplot, i_time, phases.dispersion + i*360., psym=3, $ color=colour[0] + colour[1]*(abs(i+1) mod 2)*20. ; PLOT LINES AT MULTIPLES OF 360 DEGREES, I.E. THE WRAPPING POINTS ;------------------------------------------------------------------------------- if !d.name eq 'PS' then colour = 230 else colour = 50 for i=-3,2 do oplot, [0,1000], i*360.+[180,180], color=colour ; PLOT REFERENCE PHASES IF ANY WERE PROVIDED ;------------------------------------------------------------------------------- loadct, 39, /silent if size(refphi, /type) eq 8 then begin tmpvar = (refphi.time[0]-refphi[0].time[0]) * 24.*3600. oplot, tmpvar, refphi.dispersion, color=254, thick=2 endif ; CHECK IF FLAGGING INFORMATION IS PRESENT OR NOT ;------------------------------------------------------------------------------- if keyword_set(noflag) then begin ; IF THERE IS NO FLAGGING INFO, SIMPLY PLOT ALL PHASES AT ONCE ;----------------------------------------------------------------------- oplot, i_time, unwrap, color= 50, psym=pltsym, symsiz=symsiz endif else begin ; IF THERE IS FLAGGING INFORMATION, FIRST PLOT ALL PHASES ;----------------------------------------------------------------------- oplot, i_time, unwrap, color=100, psym=pltsym, symsize=symsiz ; DETERMINE ALL ELEMENTS WHICH ARE FLAGGED AS GOOD DATA ;----------------------------------------------------------------------- flagok = midi_flagok(flgdat, phases) ; GET BOTTOM AND TOP INDICES FOR THE GOOD DATA SECTIONS ;----------------------------------------------------------------------- flagok = [flagok, 0B] idxbot = where((flagok - shift(flagok, 1)) eq 1B); > 0 < (number-1) idxtop = where((flagok - shift(flagok, -1)) eq 1B); > 0 < (number-1) ; PLOT THE PHASES FOR ALL GOOD DATA SECTIONS ;----------------------------------------------------------------------- for i=0L,n_elements(idxbot)-1 do begin oplot, i_time[idxbot[i]:idxtop[i]], $ unwrap[idxbot[i]:idxtop[i]], color=50, $ psym=pltsym, symsize=symsiz endfor endelse ; SUPPRESS LABELS ON LOWER X AXIS (USED FOR COMBINATION WITH PHASE PLOT) ;------------------------------------------------------------------------------- if keyword_set(noxlab) then begin xtitle = '' xtickn = replicate(' ', 29) endif else begin xtitle = 'frame number' xtickn = '' endelse ; OVER PLOT THE AXIS, TITLES AND LABLES ;------------------------------------------------------------------------------- plot, xrange, yrange, yrange=yrange, xtitle='time t [sec]', $ ytitle='water vapour phase [deg]', xstyle=9, ystyle=1, $ position=pltpos, /device, /nodata, /noerase axis, xaxis=1, xrange=interpol(findgen(n_elements(i_time)), i_time, $ xrange), xstyle=1, xtickname=xtickn, xtitle=xtitle if keyword_set(psfile) then tmpvar = 1.09*pltpos[3] else tmpvar = 1.08*pltpos[3] xyouts, 0.5*total(pltpos[[0,2]]), tmpvar, titles, align=0.5, $ charsize=1.20*!p.charsize, /device ; CLOSE THE PLOT AND RESET THE PLOTTING ENVIRONMENT ;------------------------------------------------------------------------------- if size(pltset, /type) eq 8 then setplot, pltset, /close END ;+ ; NAME: ; MIDI_PLTOPD ; ; PURPOSE: ; Plot the Optical Path Difference for a MIDI fringe track. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; midi_pltopd, filnam [, xrange=xrange] [, yrange=yrange] ; [, noflag=noflag] [, titles=titles] ; [, symbol=symbol] [, pltimg=pltimg] ; [, gsmoot=gsmoot] [, asmoot=asmoot] ; [, minmax=minmax] [, invers=invers] ; [, phases=phases] [, noinst=noinst] ; [, nogdel=nogdel] [, pltpos=pltpos] ; [, noxlab=noxlab] [, psfile=psfile] ; ; REQUIRED INPUTS: ; filnam EWS tag or file name of the MIDI group delay data ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; xrange range of the x axis (time in seconds) ; yrange range of the y axis (delay in mm) ; noflag do not mark the flagged frames by a different colour ; titles string containing a custom title for the plot ; symbol mark the individual measurements with small points ; pltimg plot an image of the Fourier transform data ; gsmoot width for complex smoothing of the Fourier transform data ; asmoot width for amplitude smoothing of the Fourier transform data ; minmax minimum and maximum range of data for colour scaling ; invers invert the colour scaling of the Fourier transform data ; phases plot the phase instead of the amplitude of the fringes ; noinst do not plot the instrumental OPD ; nogdel do not plot the group delay determined by EWS ; refdel reference delay data to be plotted in red ; getdel determine the group delays and plot them ; maxopd maximum delay for searching the group delay ; pltpos manually specified device coordinates for the plot position ; noxlab suppress labeling of lower x axis ; flip reverse the image in y direction (for older EWS versions) ; psfile name of the PS file for PS output ; ; OUTPUTS: ; getdel group delays determined using this routine ; a plot of the Optical Path Difference ; ; DESCRIPTION AND EXAMPLE: ; This procedure loads the Group Delay (from *.groupdelay.fits), the ; instrumental Optical Path Difference (OPD) and the flagging infor- ; mation (from *.flag.fits) for an EWS tag specified in FILNAM in ; order to produce a plot of the three quantities. If a full filename ; is specified instead of an EWS tag, e.g. 'tag.powerdelay.fits', this ; file is used instead of tag.groupdelay.fits. The group delay of ; frames flagged as "good" is overplotted in a different colour over ; the continuous group delay unless the keyword NOFLAG is set. The ; plot is automatically scaled to contain the entire fringe track and ; the entire OPD range except for the offset at the beginning of the ; track. Custom ranges for the x and y axis can be specified with the ; XRANGE and YRANGE keywords. If the keyword PLTIMG is set, also the ; Fourier transform data is plotted in the background. Reference delay ; data can be provided using the keyword REFDEL, which must be an array ; of structures as obtained when using oirGetDelay(). The procedure can ; also be used to obtain estimates of the group delay and of the fringe ; amplitude by providing a named variable to the keyword GETDEL. First ; the brightest pixel within MAXOPD of the centre of the Fourier ; transform is identified. To obtain a more accurate estimate of the ; group delay the weighted moment of a few pixels around this value ; are calculated. An estimate of the fringe amplitude is calculated ; as the sum of the nearest 9 pixels to the group delay estimate. ; In its simplest form, the procedure is called by: ; midi_pltopd, 'tag' ; ; CALLED BY: ; midi_pltres ; ; CALLING: ; midi_flagok ; setplot ; ; MODIFICATION HISTORY: ; 2008-10-14 Adapted from an earlier version by Konrad R. Tristram ; 2009-05-12 K. Tristram: Added use of xstyle and ystyle for the plot. ; 2009-05-18 K. Tristram: Changed automatic title generation. ; 2009-06-04 K. Tristram: Added patch for YRANGE due to missing labels. ; 2010-04-14 K. Tristram: Improved automatic calculation of YRANGE. ; 2010-08-10 K. Tristram: Add support for irregular groupdelay. ; 2010-10-14 K. Tristram: Fixed problem with flagging ranges. ; 2011-06-16 K. Tristram: Major rewrite using MIDI_FLAGOK and option to ; plot the group delay data. ; 2011-10-21 K. Tristram: Fixed error when plotting images with gaps. ; 2011-11-08 K. Tristram: Make GETDEL same units as other delays. ; 2011-11-09 K. Tristram: Make smoothing similar to that of oirPowerDelay. ; 2011-12-07 K. Tristram: Fixed feedback in YRANGE from 2009-06-04 patch. ; 2012-01-18 K. Tristram: Correct error in zeropoint of group delay times. ; 2012-04-03 K. Tristram: Changed counters in loops variables to LONG. ; 2012-06-04 K. Tristram: Added keywords PLTPOS and NOXLAB. ; 2012-07-19 K. Tristram: Constrain range of group delay estimates. ; 2012-08-21 K. Tristram: Added keyword FLIP and calculation of YRANGE. ; 2012-09-17 K. Tristram: Fixed error when no good frames were found. ; PRO MIDI_PLTOPD, FILNAM, XRANGE=XRANGE, YRANGE=YRANGE, NOFLAG=NOFLAG, $ TITLES=TITLES, SYMBOL=SYMBOL, PLTIMG=PLTIMG, GSMOOTH=GSMOOT, $ ASMOOTH=ASMOOT, MINMAX=MINMAX, INVERS=INVERS, PHASES=PHASES, $ NOINST=NOINST, NOGDEL=NOGDEL, REFDEL=REFDEL, GETDEL=GETDEL, $ MAXOPD=MAXOPD, PLTPOS=PLTPOS, NOXLAB=NOXLAB, FLIP=FLIP, $ PSFILE=PSFILE ; INITIALISE A FEW VARIABLES ;------------------------------------------------------------------------------- if not keyword_set(gsmoot) then gsmoot = 0 if not keyword_set(asmoot) then asmoot = 0 if not keyword_set(maxopd) then maxopd = 150 expone = 1.0 ; exponent for colour scaling of the Fourier transform data delwid = 20.0 ; width for determining the subpixel peak position debug = 0B ; default of debugging parameter: no debugging plots ; DETERMINE IF ONLY EWSTAG WAS GIVEN OR PERHAPS A FULL FILENAME ;------------------------------------------------------------------------------- if strmatch(strmid(filnam, 4, 5, /reverse), '.fits') then begin tmpvar = strsplit(filnam, '.', /extract) ewstag = tmpvar[0] delfil = '.'+strjoin(tmpvar[1:*], '.') endif else begin ewstag = filnam delfil = '.groupdelay.fits' ; standard file for the group delay data ; delfil = '.fourier.fits' ; alternative file for the group delay data endelse ; LOAD THE GROUPDELAY AND FOURIER TRANSFORM DATA ;------------------------------------------------------------------------------- gdelay = oirgetdelay(ewstag+delfil) message, 'Loaded group delay data from '+ewstag+delfil+'.', /continue if file_test(ewstag+delfil) then fftdat=oirgetdata(ewstag+delfil, ierr=ierror) $ else ierror=1B if ierror then begin fftdat = oirgetdata(ewstag+'.fringes.fits', ierr=ierror) pltimg = 0B message, 'No FFT data found in '+ewstag+delfil+'.', /continue message, 'Loaded instrumental delay only from '+ewstag + $ '.fringes.fits instead.', /continue endif else message, 'Loaded instrumental delay and FFT data from ' + $ ewstag+delfil+'.', /continue ; REPLACE FRAMES LACKING TIME AND OPD INFORMATION BY PREVIOUS FRAME ;------------------------------------------------------------------------------- tmpvar = where(fftdat.time[0] eq 0.0) if tmpvar[0] ge 0 then fftdat[tmpvar] = fftdat[(tmpvar - 1) > 0] idelay = (fftdat.opd[1]-fftdat.localopd[0])*1e3 ; LOAD THE FLAGGING INFORMATION ;------------------------------------------------------------------------------- if (~ keyword_set(noflag)) then if file_test(ewstag+'.flag.fits') then begin tmpvar = obj_new('fitstable', ewstag + '.flag.fits', extname='FLAG') flgdat = tmpvar->readrows() obj_destroy, tmpvar message, 'Loaded flagging data from '+ewstag+'.flag.fits', /cont endif else begin message, 'No flagging data found.', /continue noflag = 1B endelse ; IF POINTS ARE DESIRED THEN CREATE A POINT AS A USER SYMBOL ;------------------------------------------------------------------------------- if keyword_set(symbol) then begin usersym, cos(findgen(17) * (!pi*2/16.)), $ sin(findgen(17) * (!pi*2/16.)), /FILL pltsym = -8 endif else pltsym = 0 ; SET THE SYMBOL SIZE DEPENDING ON IF IT IS A PS PLOT OR NOT ;------------------------------------------------------------------------------- if keyword_set(psfile) then symsiz=0.30 else symsiz=0.60 ; MAKE A TITLE FOR THE PLOT IF NONE HAS BEEN PROVIDED ;------------------------------------------------------------------------------- if ~ keyword_set(titles) then begin titles = strmid(oirGetKey(ewstag+delfil, 'DATE-OBS'), 0, 19) + $ ' - ' + oirGetKey(ewstag+delfil, 'DET NRTS MODE') + ': ' + $ oirGetKey(ewstag+delfil, 'OBS TARG NAME') endif ; CALCULATE THE TIMES AND DETERMINE THE XRANGE IN TIME AND IN COLUMNS ;------------------------------------------------------------------------------- i_time = reform(fftdat.time-fftdat[0].time[0])*24*3600 i_indx = findgen(n_elements(i_time)) g_time = reform(gdelay.time-fftdat[0].time[0])*24*3600 if (n_elements(xrange) eq 2) then begin x_pixr = interpol(i_indx, i_time, xrange) endif else begin x_pixr = [-0.5, n_elements(i_time)-0.5] xrange = interpol(i_time, i_indx, x_pixr) endelse ; DETERMINE YRANGE ;------------------------------------------------------------------------------- if (n_elements(yrange) ne 2) then begin ; Don't use the jump of 5 mm at the start: tmpvar = idelay[where(idelay gt (idelay[0] -2.5))] tmpvar = [tmpvar, reform(gdelay.delay)/1000.] tmpvar = tmpvar[sort(tmpvar)] yrange = tmpvar[[0.01,0.99]*n_elements(tmpvar)] yrange += (yrange[1]-yrange[0]) * [-0.1, +0.1] yrangc = yrange ; old code: ;yrange = [min([idelay[700:*]*1e6, reform(gdelay[700:*].delay)]), $ ; max([idelay[700:*]*1e6, reform(gdelay[700:*].delay)])] ystyle = 0 endif else begin yrangc = yrange + [-0.01,+0.01]*(yrange[1]-yrange[0]) ystyle = 1 endelse ; SET THE PLOTTING ENVIRONMENT IF NOT PROVIDED AND LOAD THE COLOUR TABLE ;------------------------------------------------------------------------------- if n_elements(pltpos) ne 4 then begin setplot, pltset, pltpos, /longer, psfile=psfile pltpos *= [1.0, 1.0, 1.0, 0.94] endif loadct, 39, /silent ; CREATE THE THE COORDINATE SYSTEM OF THE PLOT ;------------------------------------------------------------------------------- plot, g_time, gdelay.delay, xrange=xrange, xstyle=5, yrange=yrangc, $ ystyle=4+ystyle, position=pltpos, /device, /nodata, color=220, /noerase ; CHECK IF THE FOURIER TRANSFORM IS TO BE PLOTTED AS AN IMAGE ;------------------------------------------------------------------------------- if keyword_set(pltimg) then begin ; CONVERT PSEUDOCOMPLEX DATA TO IDL COMPLEX DATA ;----------------------------------------------------------------------- if n_elements(fftdat[0].data1) eq 512 then begin pltarr = complex(transpose(sqrt(fftdat.data1 > 0)), 0.0) endif else pltarr = pseudocomplex(fftdat.data1) pltarr[0:1,*] = 0.0 ; APPLY COMPLEX SMOOTHING TO THE DATA ARRAY ;----------------------------------------------------------------------- ; This is slightly different to 'complexGaussSmooth' in oirPowerDelay if gsmoot ge 0.0 then pltarr = csmooth2(pltarr, gsmoot) ; GET AMPLITUDE OR PHASE OF THE COMPLEX NUMBERS ;----------------------------------------------------------------------- if keyword_set(phases) then pltarr = atan(pltarr, /phase) $ else pltarr = abs(pltarr) ; SMOOTH THE DATA IN TIME DIRECTION ;----------------------------------------------------------------------- if asmoot ne 0.0 then if asmoot gt 0.0 then begin ; For boxcar smoothing (as used in oirPowerDelay using 'minrtsBoxSmooth'), ; use the standard IDL smoothing routine: tmpvar = asmoot * (2. * sqrt(2.*alog(2.))) pltarr = smooth(pltarr, [asmoot, 0]) endif else begin ; For Gaussian smoothing (as used implemented in oirPowerDelay using ; 'minrtsBoxSmooth') use EWS routine with a very small y width: tmpvar = asmoot; / (2. * sqrt(2.*alog(2.))) pltarr = fsmooth2(pltarr, tmpvar, 0.001) endelse ; GET THE OPD STEPPING OF THE FOURIER TRANSFORM DATA ;----------------------------------------------------------------------- tmpvar = obj_new('FITSTABLE', ewstag+delfil, extname='IMAGING_DATA') header = tmpvar->head() opdstp = 1000.*header->getpar('OPD2') opdoff = 1000.*header->getpar('OPD0') obj_destroy, tmpvar obj_destroy, header ; Need the following for older snapshots of EWS and .PowerDelay.fits: ; if opdstp eq 0.0 then begin ; tmpvar = obj_new('FITSTABLE', ewstag+'.groupdelay.fits', extname='IMAGING_DATA') ; header = tmpvar->head() ; opdstp = 1000.*header->getpar('OPD2') ; opdoff = 1000.*header->getpar('OPD0') ; obj_destroy, tmpvar ; obj_destroy, header ; endif ; REVERSE THE ARRAY IF OPD STEPPING IS NEGATIVE (WHICH IT IS NORMALLY) ;----------------------------------------------------------------------- if opdstp lt 0.0 then begin opdstp = - opdstp pltarr = reverse(pltarr, 2) endif opdoff -= (n_elements(pltarr[0,*])/2.-1) * opdstp ; HAVE TO REVERSE THE ARRAY IN Y DIRECTION FOR SOME OLDER EWS VERSIONS ;----------------------------------------------------------------------- if keyword_set(flip) then pltarr = reverse(pltarr, 2) ; DETERMINE THE Y ROWS TO BE USED FOR PLOTTING ;----------------------------------------------------------------------- y_pixr = (!y.crange - opdoff) / opdstp y_pixr = fix(y_pixr + [1.0, 0.0]) > 0 < (n_elements(pltarr[0,*])-1) ; FIND THE PEAK OF THE DELAY FUNCTION ;----------------------------------------------------------------------- if arg_present(getdel) and ~ keyword_set(phase) then begin ; PREPARE THE STRUCTURE HOLDING THE NEW DELAY ESTIMATES ;--------------------------------------------------------------- getdel = replicate(gdelay[0], n_elements(i_time)) getdel.time = fftdat.time ; DETERMINE THE RANGE OF PIXELS CORRESPONDING TO MAXOPD ;--------------------------------------------------------------- pixran = 0.002 * maxopd / opdstp * [-1, +1] pixran = fix(pixran + n_elements(pltarr[0,*])/2. + 0.5) pixran = pixran > 0 < (n_elements(pltarr[0,*])-1) ; LOOP OVER THE INDIVIDUAL FRAMES ;--------------------------------------------------------------- for i=0L,n_elements(i_time)-1 do begin ; IF IN DEBUGGING MODE THEN MAKE A PLOT ;------------------------------------------------------- ; if i eq 615 then debug = 1B if debug then begin xy_tmp = [!x, !y] tmpvar = indgen(pixran[1]-pixran[0]+1)+pixran[0] plot, tmpvar, pltarr[i, tmpvar], xr=[220,300], $ psym=-8, xtitl='row', ytitl='amplitude', $ title=string(i, format='("frame:", I6)') endif ; FIND THE POSITION OF THE BRIGHTEST ELEMENT ;------------------------------------------------------- tmpvar = max(pltarr[i,pixran[0]:pixran[1]], maxidx) maxidx += pixran[0] if debug then plots, maxidx, tmpvar, psym=6, $ color=250, thick=3 ; BIAS OPD POSITION TOWARDS PREVIOUSLY FOUND OPD ;------------------------------------------------------- if i gt 0 then if tmpvar lt (2*getdel[i-1].amplitude) $ then maxidx = fix((getdel[i-1].delay/1000. - $ opdoff) / opdstp) if debug then plots, maxidx, pltarr[i,maxidx], psym=7, $ color=250, thick=3, symsize=2 ; METHOD 1 TO GET A SUB PIXEL PEAK ESTIMATE ;------------------------------------------------------- ; Quite fast, but has jumps because of moving subarray tmpvar = findgen(delwid+1) - delwid/2. + maxidx[0] if 0 then begin if debug then oplot, tmpvar, pltarr[i,tmpvar]-0E7, $ thick=2, color=200 maxidx = min(pltarr[i,tmpvar]) maxidx = total((pltarr[i,tmpvar]-maxidx) * tmpvar) / $ total((pltarr[i,tmpvar]-maxidx)) if debug then oplot, maxidx*[1,1], [0E00,1E10], $ color=250 getdel[i].delay = 1000. * (opdoff + opdstp * maxidx) ; GET AN ESTIMATE OF THE FRINGE AMPLITUDE ;------------------------------------------------------- tmpvar = fix(maxidx + [-4.0, +4.0]) getdel[i].amplitude = $ total(pltarr[i, tmpvar[0]:tmpvar[1]]) ; METHOD 2 TO GET A SUB PIXEL PEAK ESTIMATE ;------------------------------------------------------- ; Slower, because doing Gaussian fit for every frame endif else begin gaufit = gaussfit(tmpvar, pltarr[i,tmpvar], gaupar, $ estima=[pltarr[i, maxidx], maxidx, 2., 0.], $ nterms=4) gaupar[1] = gaupar[1] > pixran[0] < pixran[1] getdel[i].delay = 1000. * (opdoff + opdstp * gaupar[1]) getdel[i].amplitude = gaupar[0] if debug then oplot, tmpvar, gaufit, thick=2, color=200 endelse ; IF IN DEBUGGING MODE THEN CLOSE PLOT ;------------------------------------------------------- if debug then begin debug = '' read, debug if debug eq '' then debug = 1B else debug = 0B !x = xy_tmp[0] !y = xy_tmp[1] erase endif endfor endif ; FIND THE MINIMUM AND THE MAXIMUM OF THE DATA ARRAY ;----------------------------------------------------------------------- if n_elements(minmax) lt 2 then begin tmpvar = pltarr[sort(pltarr)] minmax = tmpvar[[0.10,0.996]*n_elements(tmpvar)] endif ; APPLY THE COLOUR SCALING TO THE DATA ARRAY AND LOAD THE COLOUR TABLE ;----------------------------------------------------------------------- pltarr = (pltarr - minmax[0]) / (minmax[1] - minmax[0]) > 0 < 1 pltarr = 255. * pltarr^expone loadct, 0, /silent ; INVERT THE ARRAY RANGE IF CORRESPONDING KEYWORD IS SET ;----------------------------------------------------------------------- if keyword_set(invers) then pltarr = 255. - pltarr ; DETERMINE THE X COLUMNS TO BE USED FOR PLOTTING ;----------------------------------------------------------------------- x_pixr = fix(x_pixr + [1.0, 0.0]) > 0 < n_elements(i_time) ; CHECK IF ANYTHING OF THE ARRAY HAS TO BE PLOTTED AT ALL ;----------------------------------------------------------------------- if (x_pixr[1] gt x_pixr[0]) and (y_pixr[1] gt y_pixr[0]) then begin ; DETERMINE WHERE THERE ARE GAPS IN THE DATA ;----------------------------------------------------------------------- gappos = i_time - shift(i_time, 1) gappos = where((gappos gt (4.0*median(gappos))) and (gappos lt 3600.)) tmpvar = where((gappos gt x_pixr[0]) and (gappos lt x_pixr[1])) if tmpvar[0] ge 0 then x_pixr = [x_pixr[0], gappos[tmpvar], x_pixr[1]] ; LOOP OVER THE SUBARRAYS ;----------------------------------------------------------------------- for i=0L,n_elements(x_pixr)-2 do begin ; CUT OUT THE SUBARRAY TO BE PLOTTED FROM THE DATA ARRAY ;----------------------------------------------------------------------- subarr = pltarr[x_pixr[i]:(x_pixr[i+1]-1), y_pixr[0]:y_pixr[1]] ; DETERMINE THE DEVICE COORDINATES FOR THE SUBARRAY ;----------------------------------------------------------------------- x_pltr = interpol(i_time[x_pixr[i]:(x_pixr[i+1]-1)], $ i_indx[x_pixr[i]:(x_pixr[i+1]-1)], $ x_pixr[i:(i+1)]+[-0.5, -0.5]) x_pltr = interpol(pltpos[[0,2]], xrange, x_pltr) y_pltr = interpol(pltpos[[1,3]], !y.crange, y_pixr*opdstp+opdoff) ; RESIZE THE SUBARRAY FOR X- OR FOR PS-OUTPUT ;----------------------------------------------------------------------- tmpvar = 2000 ; maximum number of pixels in x-direction in PS files tmpvar = 1000 if !d.name eq 'X' then begin subarr = congrid(subarr, x_pltr[1]-x_pltr[0], $ y_pltr[1]-y_pltr[0]) endif else if n_elements(subarr[*,0]) gt tmpvar then begin subarr = congrid(subarr, tmpvar, n_elements(subarr[0,*])) endif ; PLOT THE SUBARRAY ;----------------------------------------------------------------------- tv, subarr, x_pltr[0], y_pltr[0], xsize=x_pltr[1]-x_pltr[0], $ ysize=y_pltr[1]-y_pltr[0] ; END OF THE LOOP OVER SUBARRAYS ;----------------------------------------------------------------------- endfor ; END OF CHECK IF ARRAY HAS TO BE PLOTTED ;----------------------------------------------------------------------- endif loadct, 39, /silent ; MARK THE PEAKS OF THE DELAY FUNCTION ;----------------------------------------------------------------------- if arg_present(getdel) and ~ keyword_set(phase) then begin oplot, i_time, getdel.delay/1000., psym=pltsym, $ symsize=symsiz, color=200 endif endif ; PLOT A REFERENCE DELAY ;------------------------------------------------------------------------------- if size(refdel, /type) eq 8 then begin tmpvar = (refdel.time[0]-refdel[0].time[0]) * 24.*3600. oplot, tmpvar, refdel.delay/1000., color=254, thick=1 endif ; OVERPLOT THE FLAGGED OPDS WITH A DIFFERENT COLOR IF NOT SUPPRESSED ;------------------------------------------------------------------------------- if keyword_set(noflag) then begin ; PLOT THE INSTRUMENTAL OPD ;----------------------------------------------------------------------- if ~ keyword_set(noinst) then $ oplot, i_time, idelay, color=150, psym=pltsym, symsize=symsiz ; PLOT THE ENTIRE GROUP DELAY IN DARK BLUE ;----------------------------------------------------------------------- if ~ keyword_set(nogdel) then $ oplot, g_time, gdelay.delay/1000., color= 50, psym=pltsym, symsiz=symsiz endif else begin ; PLOT THE INSTRUMENTAL OPD ;----------------------------------------------------------------------- if ~ keyword_set(noinst) then $ oplot, i_time, idelay, color=150, psym=pltsym, symsize=symsiz ; DETERMINE ALL ELEMENTS WHICH ARE FLAGGED AS GOOD DATA ;----------------------------------------------------------------------- flagok = midi_flagok(flgdat, fftdat) ; GET BOTTOM AND TOP INDICES FOR THE GOOD DATA SECTIONS ;----------------------------------------------------------------------- flagok = [flagok, 0B] idxbot = where((flagok - shift(flagok, 1)) eq 1B); > 0 < (number-1) idxtop = where((flagok - shift(flagok, -1)) eq 1B); > 0 < (number-1) ; PLOT THE GROUP DELAY FOR ALL GOOD DATA SECTIONS ;----------------------------------------------------------------------- loadct, 08, /silent if ~ keyword_set(noinst) then for i=0L,n_elements(idxbot)-1 do begin oplot, i_time[idxbot[i]:idxtop[i]], $ idelay[idxbot[i]:idxtop[i]], color=125, $ psym=pltsym, symsize=symsiz endfor loadct, 39, /silent ; PLOT THE ENTIRE GROUP DELAY IN A LIGHTER BLUE ;----------------------------------------------------------------------- if ~ keyword_set(nogdel) then $ oplot, g_time, gdelay.delay/1000., color=100, psym=pltsym, symsiz=symsiz ; DETERMINE ALL ELEMENTS WHICH ARE FLAGGED AS GOOD DATA ;----------------------------------------------------------------------- flagok = midi_flagok(flgdat, gdelay) ; GET BOTTOM AND TOP INDICES FOR THE GOOD DATA SECTIONS ;----------------------------------------------------------------------- flagok = [flagok, 0B] idxbot = where((flagok - shift(flagok, 1)) eq 1B); > 0 < (number-1) idxtop = where((flagok - shift(flagok, -1)) eq 1B); > 0 < (number-1) ; PLOT THE GROUP DELAY FOR ALL GOOD DATA SECTIONS ;----------------------------------------------------------------------- if (~ keyword_set(nogdel)) and (idxbot[0] ge 0) then $ for i=0L,n_elements(idxbot)-1 do begin oplot, g_time[idxbot[i]:idxtop[i]], $ gdelay[idxbot[i]:idxtop[i]].delay/1000., color= 50, $ psym=pltsym, symsize=symsiz endfor endelse ; SUPPRESS LABELS ON LOWER X AXIS (USED FOR COMBINATION WITH PHASE PLOT) ;------------------------------------------------------------------------------- if keyword_set(noxlab) then begin xtitle = '' xtickn = replicate(' ', 29) endif else begin xtitle = 'time t [sec]' xtickn = '' endelse ; OVER PLOT THE AXIS, TITLES AND LABLES ;------------------------------------------------------------------------------- plot, g_time, gdelay.delay, xrange=xrange, xstyle=9, yrange=yrangc, $ ystyle=ystyle, xtitle=xtitle, ytitle='delay d [mm]', $ xtickname=xtickn, position=pltpos, /device, /nodata, /noerase axis, xaxis=1, xrange=interpol(findgen(n_elements(i_time)), i_time, $ xrange), xstyle=1, xtitle='frame number' if keyword_set(psfile) then tmpvar = 1.09*pltpos[3] else tmpvar = 1.08*pltpos[3] xyouts, 0.5*total(pltpos[[0,2]]), tmpvar, titles, align=0.5, $ charsize=1.00*!p.charsize, /device ; CLOSE THE PLOT AND RESET THE PLOTTING ENVIRONMENT ;------------------------------------------------------------------------------- if size(pltset, /type) eq 8 then setplot, pltset, /close END ;+ ; NAME: ; MIDI_PLTOVL ; ; PURPOSE: ; Plot the beam overlap for the photometry, fringes and masks. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; midi_pltovl, dattag [, wavran] [, mskfil=mskfil] [, /normal] ; [, pltpos=pltpos] [, psfile=psfile] ; [, ierror=ierror] ; ; REQUIRED INPUTS: ; dattag data set or EWS tag of the data to be plotted ; ; OPTIONAL INPUTS: ; wavran wavelength range over which the profiles are averaged ; ; KEYWORDS: ; mskfil name of the mask file if no DATSET structure is provided ; normal normalise all the profiles ; pltpos positions of the two plot windows ; psfile name of the PS file for PS output ; ; OUTPUTS: ; plot of the overlap of the beams either on screen or as a PS file ; ; DESCRIPTION AND EXAMPLE: ; This programme reads in the photometry files (tag.Aphotometry.fits and ; tag.Bphotometry.fits) and, if present, the fringe image file ; (tag.fringeimageave.fits) for the given tag or data set. If the data ; set also contains a valid mask file name or a valid mask file name is ; specified by the MSKFIL keyword, also the mask data is read in. The ; profiles of the photometries, the correlated flux and the mask are ; then plotted as a function of the pixels in y-direction for the two ; detector windows. Note that the correlated flux profile is the same ; for both windows as the two windows have been subtracted to create ; the fringe image. If the keyword NORMAL is set, the profiles are all ; normalised to 1 at their peak. The programme is called by: ; midi_pltovl, datset, [11,13], /normal, psfile='/tmp/overlap' ; ; CALLED BY: ; midi_pltres ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2010-08-13 Written by Konrad R. Tristram ; 2011-10-21 K. Tristram: Added output of sky mask. ; 2012-02-05 K. Tristram: Fixed error when sky mask not provided. ; 2012-02-05 K. Tristram: Fixed wrong colour of photometries in window 2. ; 2012-06-06 K. Tristram: Major rewrite and addition of keyword PLTPOS. ; 2012-08-06 K. Tristram: Add output of sky bands for data set input. ; 2012-11-01 K. Tristram: Adjust PS font size for !p.font=0. ; PRO MIDI_PLTOVL, DATTAG, WAVRAN, MSKFIL=MSKFIL, NORMAL=NORMAL, PLTPOS=PLTPOS, $ PSFILE=PSFILE, IERROR=IERROR ; INITIALISE A FEW VARIABLES ;------------------------------------------------------------------------------- if n_elements(mskfil) lt 1 then mskfil = '' if n_elements(skyfil) lt 1 then skyfil = '' skypos = [-1, -1] ; CHECK IF DATTAG IS A STRUCT OR A TAG AND COPY VARIABLES ACCORDINGLY ;------------------------------------------------------------------------------- if size(dattag, /type) eq 8 then begin ewstag = dattag.ewstag mskfil = dattag.redpar.mskfil skyfil = dattag.redpar.skyfil skypos = dattag.redpar.skypos endif else if size(dattag, /type) eq 7 then ewstag = dattag $ else begin message, 'DATTAG has to be a datset or a tag name.', /continue ierror = 1B return endelse ; SET DEFAULT VALUE FOR THE WAVELENGTH RANGE TO BE AVERAGED ;------------------------------------------------------------------------------- if n_elements(wavran) ne 2 then wavran = [8.25, 9.25] ; CHECK WHICH DATA ACTUALLY EXISTS AND IS VALID ;------------------------------------------------------------------------------- filnam = ewstag+['.fringeimageave.fits', '.fringeimages.fits', $ '.Aphotometry.fits', '.Bphotometry.fits'] filidx = where(file_test(filnam)) ; Because EWS can also create corrupt files actually have to load everything ; to check if it really is ok. if filidx[0] ge 0 then for i=0,n_elements(filidx)-1 do begin tmpvar = oirgetdata(filnam[filidx[i]], ierr=ierror) if ierror then filidx[i] = -1 endfor select = where(filidx ge 0) if select[0] lt 0 then begin message, 'No data found. Returning.', /continue ierror = 1B return endif else filidx = filidx[select] ; GET THE INDEX OF THE ELEMENTS IN THE WAVELENGTH RANGE ;------------------------------------------------------------------------------- wavele = oirgetwavelength(filnam[filidx[0]]) select = where((wavele ge wavran[0]) and (wavele le wavran[1])) ; GET THE OBJECT NAME ;------------------------------------------------------------------------------- objnam = strtrim(oirGetKey(filnam[filidx[0]], 'OBS TARG NAME'), 2) ; CREATE THE STRUCTURE HOLDING THE PROFILES ;------------------------------------------------------------------------------- tmpvar = n_elements((oirGetData(filnam[filidx[0]])).data1[0,*]) profil = replicate({n:'', d:fltarr(tmpvar), f:0B, c:0}, 3 ,2) profil.n = [['fringes, win1', 'phot A, win1', 'phot B, win1'], $ ['fringes, win2', 'phot A, win2', 'phot B, win2']] profil.c = [[30, 100, 70], [30, 100, 70]] ; TRY TO LOAD THE FRINGE IMAGE ;------------------------------------------------------------------------------- if filidx[0] le 1 then begin tmpvar = oirgetdata(filnam[filidx[0]], ierr=ierror) endif else ierror = 1B ; CHECK IF THE DATA WAS LOADED SUCCESSFULLY ;------------------------------------------------------------------------------- if ~ ierror then begin ; CHECK IF A FRINGE IMAGE EXISTS AND COPY THE DATA INTO THE STRUCTURE ;----------------------------------------------------------------------- profil[0,0].d = total(abs((pseudocomplex(tmpvar.data1))[*, select]),2)/$ n_elements(select) profil[0,0].f = 1 if total(tag_names(tmpvar) eq 'DATA2') then begin tmpvar = abs((pseudocomplex(tmpvar.data2))[*, select]) profil[0,1].d = total(tmpvar, 2)/n_elements(select) profil[0,1].f = 1 endif else profil[0,1].d = profil[0,0].d endif ; GET THE A PHOTOMETRY AND COPY THE COMPRESSED CUTS INTO THE STRUCT ;------------------------------------------------------------------------------- if total(filidx eq 2) eq 1 then begin tmpvar = oirGetData(filnam[2]) profil[1,0].d = total(tmpvar[0].data1[select, *], 1)/n_elements(select) profil[1,0].f = 1 profil[1,1].d = total(tmpvar[0].data2[select, *], 1)/n_elements(select) profil[1,1].f = 1 endif ; GET THE B PHOTOMETRY AND COPY THE COMPRESSED CUTS INTO THE STRUCT ;------------------------------------------------------------------------------- if total(filidx eq 3) eq 1 then begin tmpvar = oirGetData(filnam[3]) profil[2,0].d = total(tmpvar[0].data1[select, *], 1)/n_elements(select) profil[2,0].f = 1 profil[2,1].d = total(tmpvar[0].data2[select, *], 1)/n_elements(select) profil[2,1].f = 1 endif ; NORMALISE THE DATA IF DESIRED ;------------------------------------------------------------------------------- if keyword_set(normal) then for i=0,5 do $ if profil[i].f then profil[i].d /= max(profil[i].d[5:25]) ; CALCULATE THE X-VALUES AND THE Y-RANGE OF THE PLOT ;------------------------------------------------------------------------------- xvalue = findgen(n_elements(profil[0,0].d)) yrange = [min(profil.d[5:25]), 1.01*max(profil.d[5:25])] ; SET PARAMETERS OF THE Y-AXIS ACCORDING TO THE MODE USED ;------------------------------------------------------------------------------- if keyword_set(normal) then begin yrange = fix(10.0 * yrange - [1.0, -2.0]) / 10. > (-1.0) ystyle = 1 yunits = 'normalised' endif else begin ystyle = 0 yunits = 'ADU' endelse ; CHECK IF A THE OBJECT MASK EXISTS AND COPY THE DATA INTO A VARIABLE ;------------------------------------------------------------------------------- if file_test(mskfil) then begin tmpvar = oirgetdata(mskfil) mskdat = fltarr(n_elements(tmpvar[0].data1[0,*]), 2) mskdat[*,0] = total(tmpvar[0].data1[select, *], 1)/n_elements(select) mskdat[*,1] = total(tmpvar[0].data2[select, *], 1)/n_elements(select) if keyword_set(normal) then begin for i=0,1 do mskdat[*,i] /= max(mskdat[5:25,i]) endif else mskdat *= yrange[1] endif ; CHECK IF A THE SKY MASK EXISTS AND COPY THE DATA INTO A VARIABLE ;------------------------------------------------------------------------------- if file_test(skyfil) then begin tmpvar = oirgetdata(skyfil) skydat = fltarr(n_elements(tmpvar[0].data1[0,*]), 2) skydat[*,0] = total(tmpvar[0].data1[select, *], 1)/n_elements(select) skydat[*,1] = total(tmpvar[0].data2[select, *], 1)/n_elements(select) if keyword_set(normal) then begin for i=0,1 do skydat[*,i] /= max(skydat[5:25,i]) endif else skydat *= yrange[1] endif ; SAVE THE OLD PLOT SETTINGS ;------------------------------------------------------------------------------- d_orig = !d.name p_orig = !p ; RESCALE THE Y AXIS IF MORE THAN 1000 COUNTS ;------------------------------------------------------------------------------- if yrange[1] ge 1000 then begin factor = 0.001 ystrin = 'k' endif else begin factor = 1.0 ystrin = '' endelse ; PREPARE THE TITLE AND LABEL OF THE Y AXIS ;------------------------------------------------------------------------------- if n_elements(pltpos) ne 8 then begin titles = 'Overlap of beams and mask for '+objnam+' ' ytitle = 'intensity I ['+ystrin+yunits+']' endif ; OPEN A NEW PLOT FILE OR WINDOW AND SET THE PLOT ENVIRONMENT ;------------------------------------------------------------------------------- if n_elements(pltpos) eq 8 then noerase=1 else begin if keyword_set(psfile) then begin set_plot, 'PS' device, file=psfile+'.ps', bits_per_pixel=8, xoffset=6.1, $ yoffset=8.3, xsize=8.8, ysize=13.2, /portrait, /color if !p.font eq 0 then !p.charsize=0.65 else !p.charsize=0.75 symsiz = 0.5 if n_elements(pltpos) ne 8 then $ pltpos = [1665, 6860, 8300, 12670, 1665, 1050, 8300, 6860] endif else begin if !version.os_family eq 'Windows' then set_plot, 'WIN' $ else set_plot, 'X' window, /free, xsize=400, ysize=600 !p.charsize=1.25 symsiz = 1.0 if n_elements(pltpos) ne 8 then $ pltpos = [75, 313, 378, 575, 75, 50, 378, 313] endelse endelse !p.multi = [0, 1, 2] ; SET THE COLOURS OF THE SHADED SKY BANDS ;------------------------------------------------------------------------------- if !d.name eq 'PS' then colour = 247 else colour = 40 ; PREPARE THE PLOT IN THE UPPER PANEL: DETECTOR WINDOW 1 ;------------------------------------------------------------------------------- plot, [0, 30], factor*yrange, xstyle=5, ystyle=4+ystyle, ymargin=[1,2], $ /nodata, position=pltpos[0:3], /noerase, /device ; MARK THE DEFAULT SKY BANDS ;------------------------------------------------------------------------------- loadct, 03, /silent if skypos[0] ge 0 then for i=0,1 do begin polyfill, skypos[i]+[-2.5, +2.5, +2.5, -2.5], $ !y.crange[[0,0,1,1]], color=colour endfor loadct, 39, /silent ; CALCULATE THE LEGEND PARAMETERS AND PLOT INFORMATION FOR THE UPPER PANEL ;------------------------------------------------------------------------------- xlegen = !x.crange[0] + [0.73, 0.78, 0.80, 0.05] * (!x.crange[1] - !x.crange[0]) ylegen = 0.08 * (!y.crange[1] - !y.crange[0]) legidx = 1.0 xyouts, xlegen[3], !y.crange[1]-(legidx+0.2)*ylegen, 'Window 1' ; LOOP OVER THE ELEMENTS IN THE PROFILE STRUCTURE ;------------------------------------------------------------------------------- for i=0,2 do if profil[i,0].f then begin ; PLOT THE PROFILE ;----------------------------------------------------------------------- oplot, xvalue, factor*profil[i,0].d, color=profil[i,0].c, thick=(2-i)>1 ; PLOT THE CORRESPONDING LEGEND ENTRY ;----------------------------------------------------------------------- oplot, xlegen[0:1], !y.crange[1]-(legidx-0.2)*ylegen*[1,1], $ color=profil[i,0].c, thick=(2-i)>1 tmpvar = (strsplit(profil[i,0].n, ',', /extract))[0] xyouts, xlegen[2], !y.crange[1]-(legidx+0.2)*ylegen, tmpvar legidx += 1.0 endif ; OVERPLOT THE MASK PROFILES FOR THE UPPER PANEL (INCLUDING LEGEND) ;------------------------------------------------------------------------------- if n_elements(mskdat) gt 0 then begin oplot, xvalue, factor*mskdat[*,0], color=250, thick=2 oplot, xlegen[0:1], !y.crange[1]-(legidx-0.2)*ylegen * [1,1], $ color=250, thick=2 xyouts, xlegen[2], !y.crange[1]-(legidx+0.2)*ylegen, 'mask' legidx += 1.0 endif if n_elements(skydat) gt 0 then oplot, xvalue, factor*skydat[*,0], color=210 ; PLOT THE AXES FOR THE UPPER PANEL ;------------------------------------------------------------------------------- plot, [0, 30], factor*yrange, xstyle=1, xtickname=replicate(' ', 30), $ ystyle=ystyle, ymargin=[1,2], ytitle=ytitle, title=titles, /nodata, $ position=pltpos[0:3], /noerase, /device ; PREPARE THE PLOT IN THE LOWER PANEL: DETECTOR WINDOW 2 ;------------------------------------------------------------------------------- plot, [0, 30], factor*yrange, xrange=xrange, xstyle=5, ymargin=[4,-1], $ ystyle=4+ystyle, /nodata, position=pltpos[4:7], /noerase, /device ; MARK THE DEFAULT SKY BANDS ;------------------------------------------------------------------------------- loadct, 03, /silent if skypos[0] ge 0 then for i=0,1 do begin polyfill, skypos[i]+[-2.5, +2.5, +2.5, -2.5], $ !y.crange[[0,0,1,1]], color=colour endfor loadct, 39, /silent ; CALCULATE THE LEGEND PARAMETERS AND PLOT THE LEGEND FOR THE LOWER PANEL ;------------------------------------------------------------------------------- xlegen = !x.crange[0] + [0.05, 0.93] * (!x.crange[1] - !x.crange[0]) xyouts, xlegen[0], !y.crange[1] - 1.2*ylegen, 'Window 2' if !p.font gt -1 then strvar = ['!9l!X', '!9m!Xm'] $ else strvar = ['!4k!X', '!4l!Xm'] xyouts, xlegen[1], !y.crange[1] - 1.2*ylegen, align=1, $ string(wavran, format='(F4.1, "'+strvar[1]+' < '+strvar[0]+ $ ' < ", F4.1, "'+strvar[1]+'")') ; LOOP OVER THE ELEMENTS IN THE PROFILE STRUCTURE ;------------------------------------------------------------------------------- for i=0,2 do if profil[i,1].f then begin ; PLOT THE PROFILE ;----------------------------------------------------------------------- oplot, xvalue, factor*profil[i,1].d, color=profil[i,1].c, thick=(2-i)>1 endif ; PLOT THE MASK PROFILES FOR THE LOWER PANEL ;------------------------------------------------------------------------------- if n_elements(mskdat) gt 0 then oplot, xvalue, factor*mskdat[*,1], color=250, $ thick=2 if n_elements(skydat) gt 0 then oplot, xvalue, factor*skydat[*,1], color=210 ; PLOT THE AXES FOR THE UPPER PANEL ;------------------------------------------------------------------------------- plot, [0, 30], factor*yrange, xrange=xrange, xstyle=1, ymargin=[4,-1], $ ystyle=ystyle, xtitle='pixels on detector in y direction', $ ytitle=ytitle, psym=8, /nodata, position=pltpos[4:7], /noerase, /device ; CLOSE THE PLOTTING DEVICE AND RESET PLOT PARAMETERS ;------------------------------------------------------------------------------- if keyword_set(psfile) then device, /close set_plot, d_orig !p = p_orig ierror = 0B END ;+ ; NAME: ; MIDI_PLTRES ; ; PURPOSE: ; Plot a page of diagnostic plots for reduced MIDI data. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; midi_pltres, dattag, psfile [, adu2jy] [, mskfil=mskfil] ; [, yrange=yrange] ; ; REQUIRED INPUTS: ; dattag data set or EWS tag of the reduced data to be plotted ; psfile name of the PS file ; ; OPTIONAL INPUTS: ; adu2jy structure containing the ADU to Jy conversion factors ; ; KEYWORDS: ; mskfil name of the mask file if no DATSET structure is provided ; yrange range of the y axis for the delay plot (in mm) ; ; OUTPUTS: ; several plots that give an overview over a set of reduced MIDI data ; ; DESCRIPTION AND EXAMPLE: ; This programme composes several diagnostic plots which can be used to ; assess the quality of MIDI data and of the data reduction process. ; For this purpose the programme calls several plotting routines from ; the MIDI package. The plots are from top to bottom and left to right: ; (1) optical path difference as a function of time ; (2) water vapour phase as a function of time ; (3) image of detector window 1 of the fringes and photometries ; (4) raw masked and total flux ; (5) raw correlated flux ; (6) raw differential phase and visibility ; (7) position of the measurement in the UV plane ; (8) beam overlap (mask, correlated flux, photometries) ; (9) conversion factors of nearby calibrators ; If the necessary files are not present then empty plots are created. ; Due to the large amount of information only PS output is possible. The ; programme is called by: ; midi_pltres, datset[0], '/tmp/pltres0' ; ; CALLED BY: ; none ; ; CALLING: ; midi_getkey ; midi_pltopd ; midi_pltphi ; midi_outdis ; midi_wavaxi ; midi_pltovl ; midi_a2jave ; ; MODIFICATION HISTORY: ; 2012-12-12 Written by Konrad R. Tristram ; PRO MIDI_PLTRES, DATTAG, PSFILE, ADU2JY, MSKFIL=MSKFIL, YRANGE=YRANGE ; CHECK IF DATTAG IS A STRUCTURE OR A STRING AND COPY VARIABLES ACCORDINGLY ;------------------------------------------------------------------------------- if size(dattag, /type) eq 8 then begin ewstag = dattag.ewstag mskfil = dattag.redpar.mskfil skyfil = dattag.redpar.skyfil commen = dattag.commen endif else if size(dattag, /type) eq 7 then begin ewstag = dattag commen = '' endif else begin message, 'DATTAG has to be a datset or a tag name.', /continue return endelse ; CHECK IF A NAME FOR THE POSTSCRIPT FILE WAS PROVIDED ;------------------------------------------------------------------------------- if ~ size(psfile, /type) eq 7 then begin message, 'A postscript file name has to be provided as the ' + $ 'second argument.', /continue stop endif else message, 'Writing output to ' + psfile + '.ps', /continue ; GET SEVERAL HEADER KEYWORDS ;------------------------------------------------------------------------------- keynam = ['OBS TARG NAME', 'MJD-OBS', 'DATE-OBS', 'COU GUID MODE', $ 'COU GUID STATUS', 'DET DIT', 'DET NRTS MODE', 'ISS AIRM END', $ 'ISS AIRM START', 'ISS AMBI FWHM END', 'ISS AMBI FWHM START', $ 'ISS AMBI TAU0 END', 'ISS AMBI TAU0 START', 'ISS CHOP ST', $ 'ISS CONF STATION1', 'ISS CONF STATION2', 'ISS IAS IRIS DIT', $ 'ISS IAS IRIS_GUID','ISS IAS SELECTED', 'OBS PROG ID'] filnam = ['.calcorr.', '.corr.', '.calvis.', '.redcal.', '.groupdelay.', $ '.ungroupdelay.', '.fringeimages.', '.fringeimageave.', $ '.calphot.', '.photometry.','.Aphotometry.', '.Bphotometry.'] tmpvar = where(file_test(ewstag+filnam+'fits') eq 1) if tmpvar[0] ge 0 then begin keywrd = midi_getkey(ewstag+filnam[tmpvar[0]]+'fits', keynam) for i=0,n_tags(keywrd)-1 do if (size(keywrd.(i), /type) eq 7) then $ keywrd.(i) = strtrim(keywrd.(i), 2) endif else begin keywrd = {obstargname:'no name', mjdobs:0.D, dateobs:'', $ couguidmode:'none', couguidstatus:'F', detdit:0.D, $ detnrtsmode:'none', issairmend:0.0, issairmstart:0.0, $ issambifwhmend:0.0, issambifwhmstart:0.0, $ issambitau0end:0.0, issambitau0start:0.0, isschopst:'F', $ isschopfreq:0.0, issconfstation1:'', issconfstation2:'', $ issiasirisdit:0.0, issiasiris_guid:'OFF', $ issiasselected:'none', obsprogid:'000.A-0000(A)'} endelse if size(dattag, /type) eq 8 then keywrd.obstargname = dattag.object ; DEFINE ADDITIONAL STRING FOR PS OUTPUT ;------------------------------------------------------------------------------- psstri = 'systemdict /setdistillerparams known { ' + string(10b) + $ '<< /AutoFilterColorImages false /ColorImageFilter ' + $ string(10b) + '/FlateEncode >> setdistillerparams} if' ; SAVE THE OLD PLOT SETTINGS IN A PLTSET STRUCTURE ;------------------------------------------------------------------------------- pltset = {d_old:!d.name, p:!p, d_new:!d.name} ; SET DEVICE AND PLOT PROPERTIES FOR THE PS OUTPUT ;------------------------------------------------------------------------------- set_plot, 'PS' device, file=psfile+'.ps', xsize=18.0, ysize=24.0, $ xoffset=1.5, yoffset=2.8, /portrait, bits_per_pixel=8, $ output=psstri, encapsulated=encaps, color=1 !p.charsize=0.5 ;!p.font = -1 ; PLOT THE DEVICE COORDINATES (FOR DEBUGGING) ;------------------------------------------------------------------------------- if 0 then plot, [0,18], [0,24], xstyle=1, ystyle=1, color=80, $ position=[0, 0, 18000, 24000], /device, /nodata, /noerase ; PLOT OUT SOME TEXT INFORMATION FROM THE HEADER ;------------------------------------------------------------------------------- xyouts, 11000, 8500, keywrd.obstargname, $ charsize=2.5*!p.charsize, charthick=3, /device xyouts, 11000, 8150, keywrd.detnrtsmode, $ charsize=1.5*!p.charsize, /device xyouts, 11000, 7800, strmid(ewstag, 0, 39), charsize=1.5*!p.charsize, /device xyouts, 11300, 500, charsize=1.25*!p.charsize, orientation=90, /device, $ 'Date: ' + strmid(keywrd.dateobs, 0, 19) + '; DETDIT =' + $ string(1000.*keywrd.detdit, format='(F5.1)') + ' ms' xyouts, 11600, 500, charsize=1.25*!p.charsize, orientation=90, /device, $ 'Guiding: ' + keywrd.couguidmode + ' - ' + keywrd.couguidstatus + $ ' & ' + strtrim(keywrd.issiasselected, 2) + ' - ' + $ strtrim(keywrd.issiasiris_guid, 2) + ', DIT = ' + $ string(keywrd.issiasirisdit, format='(F5.2)') + ' sec' xyouts, 11900, 500, charsize=1.25*!p.charsize, orientation=90, /device, $ 'Airmass: ' + string(keywrd.issairmend, format='(F5.2)') + $ '; seeing: ' + string(keywrd.issambifwhmend, format='(F5.2)') + $ ' arcsec; tau0: ' + string(1000.*keywrd.issambitau0end, $ format='(F5.2)') + ' ms' xyouts, 12200, 500, charsize=1.25*!p.charsize, orientation=90, /device, $ 'Comment: '+strmid(commen, 0, 40) ; SET THE LENGTH OF THE TICK MARKS ;------------------------------------------------------------------------------- !x.ticklen = 0.040 !y.ticklen = 0.008 ; CHECK IF A ***GROUPDELAY*** FILE EXISTS ;------------------------------------------------------------------------------- if file_test(ewstag+'.groupdelay.fits') then begin ; TRY TO ALSO GET THE POWERDELAY ;----------------------------------------------------------------------- if file_test(ewstag+'.powerdelay.fits') then begin tmpvar = oirgetdelay(ewstag + '.powerdelay.fits') if size(tmpvar, /type) eq 8 then begin select = where(tmpvar.time ne 0.) if select[0] ge 0 then tmpvar = tmpvar[select] endif else tmpvar = 0 endif else tmpvar = 0 ; MAKE THE GROUPDELAY PLOT ;----------------------------------------------------------------------- midi_pltopd, ewstag, /pltimg, gsmooth=3, asmooth=-10, titles=' ', $ pltpos=[1100, 20000, 17800, 23500], refdel=tmpvar, $ /invers, /noxlab, yrange=yrange endif else begin ; MAKE AN EMPTY PLOT ;----------------------------------------------------------------------- plot, [0.0,167.3], [+0.0, +1.0], xstyle=1, ystyle=1, /nodata, $ xtickname=replicate(' ', 29), $ positi=[1100, 20000, 17800, 23500], /device, /noerase xyouts, 9450, 21500, 'No group delay data found.', $ charsize=2.0*!p.charsize, align=0.5, /device endelse ; CHECK IF AN ***UNGROUPDELAY*** FILE EXISTS ;------------------------------------------------------------------------------- if file_test(ewstag+'.ungroupdelay.fits') then begin ; MAKE THE PHASE PLOT ;----------------------------------------------------------------------- midi_pltphi, ewstag, smtval=10, yrange=[-920, +920], titles=' ', $ pltpos=[1100, 16500, 17800, 20000], /noxlab endif else begin ; MAKE AN EMPTY PLOT ;----------------------------------------------------------------------- plot, [0.0,167.3], [-920, +920], xstyle=1, ystyle=1, /nodata, $ xtitle='time t [sec]', $ positi=[1100, 16500, 17800, 20000], /device, /noerase xyouts, 9450, 18000, 'No water vapour phase data found.', $ charsize=2.0*!p.charsize, align=0.5, /device endelse ; PRINT OUT THE NUMBER OF GOOD FRAMES ;------------------------------------------------------------------------------- goodst = total(midi_flagok(ewstag)) xyouts, 1300, 16800, string(goodst, format='(I6)'), $ charsize=2.0*!p.charsize, align=0.0, /device ; SET A FEW ADDITIONAL VARIABLES NEEDED FOR PLOTTING THE DETECTOR WINDOWS ;------------------------------------------------------------------------------- pscale = 0.086 offset = 18 arrsiz = [171, 41] ; DETERMINE THE TRACE OF THE MASK ;------------------------------------------------------------------------------- mskpos = replicate(-1, arrsiz[0]) if (size(mskfil, /type) eq 7) then if file_test(mskfil) then begin tmpdat = (oirgetdata(mskfil)).data1 mskpos = make_array(n_elements(tmpdat[*,0]), value=1)#findgen(41)*tmpdat mskpos = total(mskpos, 2)/total(tmpdat, 2) endif ; TRY TO LOAD THE ***FRINGE IMAGE*** ;------------------------------------------------------------------------------- if file_test(ewstag+'.fringeimages.fits') then begin detwin = oirgetdata(ewstag+'.fringeimages.fits', ierr=ierror) endif else ierror = 1B ; CHECK IF THE DATA WAS LOADED SUCCESSFULLY ;------------------------------------------------------------------------------- if ~ ierror then begin ; GET THE DATA OF THE DETECTOR WINDOW ;----------------------------------------------------------------------- detwin = abs(transpose(pseudocomplex(detwin.data1))) ; FIND THE MINIMUM AND THE MAXIMUM OF THE DATA ;----------------------------------------------------------------------- srtdat = sort(detwin) xrange = detwin[[srtdat[100], srtdat[n_elements(detwin)-20]]] ; RESCALE THE DATA BY THE SQUARE ROOT ;----------------------------------------------------------------------- detwin = (detwin - xrange[0]) / (xrange[1] - xrange[0]) > 0 < 1 detwin = 255 - sqrt(detwin) * 255 ; PLOT THE DATA OF THE WINDOW ;----------------------------------------------------------------------- midi_outdis, detwin, [1100, 13910, 7940, 15550], offset, $ ctable=3, disper='prism', /tlabel, /llabel, /rlabel ; OVERPLOT THE TRACE OF THE MASK ;----------------------------------------------------------------------- oplot, findgen(arrsiz[0]), mskpos, color=100, linestyle=2 endif else begin ; MAKE AN EMPTY PLOT ;----------------------------------------------------------------------- plot, [0, arrsiz[0]]-0.5, [0, arrsiz[1]]-0.5, ystyle=9, xstyle=9, $ xtickname=replicate(' ', 29), xticklen=0.08, yticklen=0.02, $ positi=[1100, 13910, 7940, 15550], /device, /noerase, /nodata midi_wavaxi, 'prism', xtickl=0.08, /pltlab yrange = ([-0.5, arrsiz[1]-0.5] - offset) * pscale axis, yaxis=1, yrange=yrange, ystyle=1, yticklen=0.02, $ ytickin=1, yminor=5 xyouts, 4520, 14650, 'No fringe image found.', $ charsize=2.0*!p.charsize, align=0.5, /device endelse xyouts, 1300, 14110, 'fringes', charsize=1.5*!p.charsize, /device ; CHECK IF A ***PHOTOMETRY A*** FILE EXISTS AND GET THE DATA ;------------------------------------------------------------------------------- if file_test(ewstag+'.Aphotometry.fits') then begin tmpvar = oirgetdata(ewstag+'.Aphotometry.fits', ierr=ierror) if not ierror then begin detwin = tmpvar[0].data1 photom = [1, 0] endif else photom = [0, 0] endif else photom = [0, 0] ; CHECK IF A ***PHOTOMETRY B*** FILE EXISTS AND GET THE DATA ;------------------------------------------------------------------------------- if file_test(ewstag+'.Bphotometry.fits') then begin tmpvar = oirgetdata(ewstag+'.Bphotometry.fits', ierr=ierror) if not ierror then begin detwin = [[[detwin]], [[tmpvar[0].data1]]] photom += [0, 1] endif endif ; SCALE THE ARRAY FOR PLOTTING ;------------------------------------------------------------------------------- if total(photom) gt 0 then begin ; FIND THE MINIMUM AND THE MAXIMUM OF THE DATA ;----------------------------------------------------------------------- srtdat = sort(detwin) xrange = detwin[[srtdat[150], srtdat[n_elements(detwin)-20]]] ; RESCALE THE DATA BY THE SQUARE ROOT ;----------------------------------------------------------------------- detwin = (detwin - xrange[0]) / (xrange[1] - xrange[0]) > 0 < 1 detwin = 255 - sqrt(detwin) * 255 endif ; MAKE A PLOT FOR ***PHOTOMETRY A*** ;------------------------------------------------------------------------------- if photom[0] eq 1 then begin ; PLOT THE DATA OF THE WINDOW ;----------------------------------------------------------------------- tmpvar = 0 midi_outdis, detwin[*,*,tmpvar], [1100, 12180, 7940, 13820], offset, $ ctable=3, disper='prism', /llabel, /rlabel ; OVERPLOT THE TRACE OF THE MASK ;----------------------------------------------------------------------- oplot, findgen(arrsiz[0]), mskpos, color=100, linestyle=2 endif else begin ; MAKE AN EMPTY PLOT ;----------------------------------------------------------------------- plot, [0, arrsiz[0]]-0.5, [0, arrsiz[1]]-0.5, ystyle=9, xstyle=9, $ xtickname=replicate(' ', 29), xticklen=0.08, yticklen=0.02, $ positi=[1100, 12180, 7940, 13820], /device, /noerase, /nodata midi_wavaxi, 'prism', xtickl=0.08 yrange = ([-0.5, arrsiz[1]-0.5] - offset) * pscale axis, yaxis=1, yrange=yrange, ystyle=1, yticklen=0.02, $ ytickin=1, yminor=5 xyouts, 4520, 12900, 'No photometry A found.', $ charsize=2.0*!p.charsize, align=0.5, /device endelse xyouts, 1300, 12380, 'phot A', charsize=1.5*!p.charsize, /device ; MAKE A PLOT FOR ***PHOTOMETRY B*** ;------------------------------------------------------------------------------- if photom[1] eq 1 then begin ; PLOT THE DATA OF THE WINDOW ;----------------------------------------------------------------------- tmpvar = photom[0] midi_outdis, detwin[*,*,tmpvar], [1100, 10450, 7940, 12090], 15, $ ctable=3, disper='prism', /llabel, /rlabel, /blabel ; OVERPLOT THE TRACE OF THE MASK ;----------------------------------------------------------------------- oplot, findgen(arrsiz[0]), mskpos, color=100, linestyle=2 endif else begin ; MAKE AN EMPTY PLOT ;----------------------------------------------------------------------- plot, [0, arrsiz[0]]-0.5, [0, arrsiz[1]]-0.5, ystyle=9, xstyle=9, $ xticklen=0.08, yticklen=0.02, $ positi=[1100, 10450, 7940, 12090], /device, /noerase, /nodata midi_wavaxi, 'prism', xtickl=0.08 yrange = ([-0.5, arrsiz[1]-0.5] - offset) * pscale axis, yaxis=1, yrange=yrange, ystyle=1, yticklen=0.02, $ ytickin=1, yminor=5 xyouts, 4520, 11220, 'No photometry B found.', $ charsize=2.0*!p.charsize, align=0.5, /device endelse xyouts, 1300, 10650, 'phot B', charsize=1.5*!p.charsize, /device ; PLOT THE AXIS LABELS FOR THE WINDOW PLOTS ;------------------------------------------------------------------------------- if !p.font gt -1 then strvar = ['!9l!X [!9m!Xm]', '!9a!X [arcsec]'] $ else strvar = ['!4k!X [!4l!Xm]', '!4a!X [arcsec]'] ;xyouts, 4520, 10000, 'pixels on detector x', align=0.5, /device xyouts, 4520, 15900, 'wavelength '+strvar[0], align=0.5, /device xyouts, 200, 13000, 'pixels on detector y', align=0.5, orient=90, /device xyouts, 8500, 13000, 'offset '+strvar[1], align=0.5, orient=90, /device ; SET A FEW ADDITIONAL VARIABLES NEEDED FOR PLOTTING THE DATA ;------------------------------------------------------------------------------- filnam = [['.photometry.fits', '.corr.fits', '.redcal.fits'], $ ['.calphot.fits', '.calcorr.fits', '.calvis.fits']] offset = 5300. * [1, 0, 1, 0] !x.ticklen = 0.04 !y.ticklen = 0.03 ; LOOP OVER RAW DATA AND CALIBRATED DATA ;------------------------------------------------------------------------------- for i=0,1 do begin ; TRY TO LOAD THE ***TOTAL FLUX*** ;------------------------------------------------------------------------------- if file_test(ewstag+filnam[0,i]) then begin tmpdat = oirgetdata(ewstag+filnam[0,i], ierr=ierror) endif else ierror = 1B ; CHECK IF THE DATA WAS LOADED SUCCESSFULLY ;------------------------------------------------------------------------------- if ~ ierror then begin ; GET THE WAVELENGTH SCALING ;----------------------------------------------------------------------- wavele = oirgetwavelength(ewstag+filnam[0,i]) ; COPY THE DATA INTO AN ARRAY ;----------------------------------------------------------------------- if i eq 0 then begin tmpdat[0].data1 = 0.5*total(tmpdat[0:1].data1, 2) tmpdat[2].data1 = sqrt(total(tmpdat[6:7].data1^2, 2)) tmpdat[1].data1 = tmpdat[5].data1 tmpdat[3].data1 = tmpdat[11].data1 endif ; DETERMINE THE SCALING ;----------------------------------------------------------------------- wavsel = where((wavele gt 7.9) and (wavele lt 13.2)) tmpvar = tmpdat[0:3].data1[wavsel] tmpvar = [tmpvar[*, 0:1] - (tmpvar[*, 2:3] 0 endif else begin ; SET DEFAULT VALUES IF NO DATA FOUND ;----------------------------------------------------------------------- ierror = 1B yrange = [0,1000] endelse ; PLOT THE ***TOTAL FLUX*** DATA ;------------------------------------------------------------------------------- plot, [7.5, 13.5], yrange, xstyle=5, ystyle=4, xtickname=replicate(' ', 29), $ positi=[1100, 6550, 5100, 9550]+i*offset, /device, /noerase, /nodata plots, [1300, 1600], [9300, 9300], color= 70, /device xyouts, 1700, 9250, 'unmasked', /device plots, [3000, 3300], [9300, 9300], color=150, /device xyouts, 3400, 9250, 'masked', /device if ~ ierror then begin oplot, wavele, tmpdat[0].data1, color=70 errplot, wavele-0.01, width=0, color=70, $ tmpdat[0].data1-tmpdat[2].data1, $ tmpdat[0].data1+tmpdat[2].data1 oplot, wavele, tmpdat[1].data1, color=150 errplot, wavele+0.01, width=0, color=150, $ tmpdat[1].data1-tmpdat[3].data1, $ tmpdat[1].data1+tmpdat[3].data1 endif else begin xyouts, 3100+i*offset[0], 8000, 'No photometry found.', $ charsize=1.8*!p.charsize, align=0.5, /device endelse plot, [7.5, 13.5], yrange, xstyle=1, ystyle=0, xtickname=replicate(' ', 29), $ positi=[1100, 6550, 5100, 9550]+i*offset, /device, /noerase, /nodata ; TRY TO LOAD THE ***CORRELATED FLUX*** ;------------------------------------------------------------------------------- if file_test(ewstag+filnam[1,i]) then begin tmpdat = oirgetvis(ewstag+filnam[1,i], wave=wavele, ierr=ierror) endif else ierror = 1B ; CHECK IF THE DATA WAS LOADED SUCCESSFULLY ;------------------------------------------------------------------------------- if ~ ierror then begin ; DETERMINE THE SCALING ;----------------------------------------------------------------------- wavsel = where((wavele gt 7.9) and (wavele lt 13.2)) tmpvar = [tmpdat.visamp[wavsel]-tmpdat.visamperr[wavsel], $ tmpdat.visamp[wavsel]+tmpdat.visamperr[wavsel]] yrange = [0, max(tmpvar)] ; yrange = [0, yr[i+2]] tmpvar = [tmpdat.visphi[wavsel]-tmpdat.visphierr[wavsel], $ tmpdat.visphi[wavsel]+tmpdat.visphierr[wavsel]] yrange = [yrange, min(tmpvar), max(tmpvar)] endif else begin ; SET DEFAULT VALUES IF NO DATA FOUND ;----------------------------------------------------------------------- ierror = 1B yrange = [0,1000, -90, 90] endelse ; PLOT THE ***CORRELATED FLUX*** DATA ;------------------------------------------------------------------------------- plot, [7.5, 13.5], yrange[0:1], xstyle=5, ystyle=4, $ xtickname=replicate(' ', 29), positi=[1100, 3550, 5100, 6550]+i*offset, $ /device, /noerase, /nodata if not ierror then begin oplot, wavele, tmpdat.visamp, color=70 errplot, wavele, tmpdat.visamp-tmpdat.visamperr, $ tmpdat.visamp+tmpdat.visamperr, width=0, color=70 endif else begin xyouts, 3100+i*offset[0], 5000, 'No corr. flux found.', $ charsize=1.8*!p.charsize, align=0.5, /device endelse plot, [7.5, 13.5], yrange[0:1], xstyle=1, ystyle=0, $ xtickname=replicate(' ', 29), $ positi=[1100, 3550, 5100, 6550]+i*offset, /device, /noerase, /nodata ; PLOT THE ***RAW PHASE*** DATA ;------------------------------------------------------------------------------- loadct, 0, /silent plot, [7.5, 13.5], yrange[2:3], xstyle=5, ystyle=4, $ xtickname=replicate(' ', 29), positi=[1100, 550, 5100, 3550]+i*offset, $ /device, /noerase, /nodata oplot, [ 7.5, 13.5], [0.0, 0.0], color=200 loadct, 39, /silent if not ierror then begin oplot, wavele, tmpdat.visphi, color=70 errplot, wavele, tmpdat.visphi-tmpdat.visphierr, $ tmpdat.visphi+tmpdat.visphierr, width=0, color=70 endif else begin xyouts, 3100+i*offset[0], 2000, 'No phase data found.', $ charsize=1.8*!p.charsize, align=0.5, /device endelse plot, [7.5, 13.5], yrange[2:3], xstyle=1, ystyle=8, $ positi=[1100, 550, 5100, 3550]+i*offset, /device, /noerase, /nodata ; TRY TO LOAD THE ***VISIBILITY*** ;------------------------------------------------------------------------------- if file_test(ewstag+filnam[2,i]) then begin tmpdat = oirgetvis(ewstag+filnam[2,i], wave=wavele, ierr=ierror) endif else ierror = 1B ; CHECK IF THE DATA WAS LOADED SUCCESSFULLY ;------------------------------------------------------------------------------- if ~ ierror then begin ; DETERMINE THE SCALING ;----------------------------------------------------------------------- wavsel = where((wavele gt 7.9) and (wavele lt 13.2)) tmpvar = [tmpdat.visamp[wavsel]-tmpdat.visamperr[wavsel], $ tmpdat.visamp[wavsel]+tmpdat.visamperr[wavsel]] yrange = [0, max(tmpvar)] < 2 endif else begin ; SET DEFAULT VALUES IF NO DATA FOUND ;----------------------------------------------------------------------- ierror = 1B yrange = [0.0,1.2] endelse ; PLOT THE ***VISIBILITY*** DATA ;------------------------------------------------------------------------------- plot, [7.5, 13.5], yrange[0:1], xstyle=5, ystyle=4, $ xtickname=replicate(' ', 29), positi=[1100, 550, 5100, 3550]+i*offset, $ /device, /noerase, /nodata if not ierror then begin oplot, wavele, tmpdat.visamp, color=150 errplot, wavele, tmpdat.visamp-tmpdat.visamperr, $ tmpdat.visamp+tmpdat.visamperr, width=0, color=150 endif axis, yaxis=1, yrange=yrange[0:1], ystyle=1, yticklen=0.02 ; END OF LOOP OVER RAW DATA AND CALIBRATED DATA ;------------------------------------------------------------------------------- endfor ; PLOT THE AXIS LABELS FOR THE RAW DATA PLOTS ;------------------------------------------------------------------------------- if !p.font gt -1 then begin strvar = ['!9l!X [!9m!Xm]', '!9f!X!Draw!N [deg]', '!9f!X!Dcal!N [deg]'] endif else begin strvar = ['!4k!X [!4l!Xm]', '!4u!X!Draw!N ['+string("232B)+']', $ '!4u!X!Dcal!N ['+string("232B)+']'] endelse xyouts, 3100, 9650, 'raw data', align=0.5, $ charsize=1.5*!p.charsize, /device xyouts, 3100+offset[0], 9650, 'calibrated data', align=0.5, $ charsize=1.5*!p.charsize, /device xyouts, 3100 , 50, 'wavelength '+strvar[0], align=0.5, /device xyouts, 3100+offset[0], 50, 'wavelength '+strvar[0], align=0.5, /device xyouts, 200 , 8050, 'raw total flux F!Dtot!N [ADU]', $ align=0.5, orientation=90, /device xyouts, 200 , 5050, 'raw corr. flux F!Dcor!N [ADU]', $ align=0.5, orientation=90, /device xyouts, 200 , 2050, 'raw phase '+strvar[1], $ align=0.5, orientation=90, /device xyouts, 450+offset[0], 8050, 'calibrated total flux F!Dtot!N [Jy]', $ align=0.5, orientation=90, /device xyouts, 450+offset[0], 5050, 'calibrated corr. flux F!Dcor!N [Jy]', $ align=0.5, orientation=90, /device xyouts, 450+offset[0], 2050, 'calibrated phase '+strvar[1], $ align=0.5, orientation=90, /device ; PLOT THE ***BEAM OVERLAP*** ;------------------------------------------------------------------------------- midi_pltovl, dattag, [11,13], ierror=ierror, /norm, $ pltpos=[13800, 12550, 17800, 15550, 13800, 9550, 17800, 12550] if ierror then begin plot, [0.0, 30.0], [-0.1, 1.2], xstyle=1, ystyle=1, $datset.ewstag xtickname=replicate(' ', 30), $ position=[13800, 12550, 17800, 15550], /device, /nodata, /noerase plot, [0.0, 30.0], [-0.1, 1.2], xstyle=1, ystyle=1, $ position=[13800, 9550, 17800, 12550], /device, /nodata, /noerase endif xyouts, 13200, 14050, 'intensity I [normalised]', align=0.5, $ orientation=90, /device xyouts, 13200, 11050, 'intensity I [normalised]', align=0.5, $ orientation=90, /device xyouts, 15800, 15650, 'beam overlap', charsize=1.5*!p.charsize, $ align=0.5, orientation=0, /device ; DEFINE A FEW VARIABLES NEED FOR PLOTTING THE TRANSFER FUNCTION ;------------------------------------------------------------------------------- wavele = [ 8.5, 12.5] ; plot the transfer function at these wavelengths colour = [ 50, 250] ; colours for the two wavelengths wavdel = 0.2 ; half width of wavelength band tagnam = 'CORFLX' ; select which tag to actually process & plot ; CHECK IF PLOTTING THE TRANSFER FUNCTIONS ;------------------------------------------------------------------------------- if (keywrd.mjdobs ne 0.0) and (size(adu2jy, /type) eq 8) then begin ; GET THE INDEX OF THE TAG THAT NEED TO BE PROCESSED ;----------------------------------------------------------------------- tagidx = where(tag_names(adu2jy) eq tagnam) ; CALCULATE THE AVERAGE AND SCATTER OF THE TRANSFER FUNCTION ;----------------------------------------------------------------------- a2jave = midi_a2jave(adu2jy, keywrd.mjdobs, 50, flgdat=flgdat, /silent) flgidx = where(tag_names(flgdat) eq tagnam) ; SELECT THE CALIBRATORS WHICH CONTAIN DATA ;----------------------------------------------------------------------- select = where((flgdat[*].(flgidx) eq 0) or (flgdat[*].(flgidx) eq 2)) a2jdat = adu2jy[select] a2jdat.mjdobs = 24. * (a2jdat.mjdobs - keywrd.mjdobs) flgdat = flgdat[select].(flgidx) ; CALCULATE THE AVERAGES IN THE WAVELENGTH BANDS ;----------------------------------------------------------------------- for i=0,n_elements(wavele)-1 do begin wavsel = where((a2jdat[0].wavele gt (wavele[i] - wavdel)) and $ (a2jdat[0].wavele lt (wavele[i] + wavdel))) for j=0,1 do begin a2jdat.(tagidx+j)[i] = total(a2jdat.(tagidx+j)[wavsel], 1)/$ n_elements(wavsel) a2jave.(tagidx+j)[i] = total(a2jave.(tagidx+j)[wavsel])/$ n_elements(wavsel) endfor endfor ; DETERMINE THE PLOT RANGES ;----------------------------------------------------------------------- xrange = [min([a2jdat.mjdobs, 0])-0.2, max([a2jdat.mjdobs, 0])+0.2] wavsel = where((a2jdat[0].wavele gt 8.0) and $ (a2jdat[0].wavele lt 13.0)) yrange = [0, max(a2jdat.(tagidx )[wavsel] + $ a2jdat.(tagidx+1)[wavsel])] ierror = 0B ; CHECK IF THE ASSOCIATED CALIBRATOR IS KNOWN ;----------------------------------------------------------------------- if size(dattag, /type) eq 8 then begin calidx = where(dattag[0].redpar.calibr eq a2jdat.identi) endif else calidx = -1 ; IF NOT PLOTTING THE TRANSFER FUNCTION PREPARE A DUMMY PLOT ;------------------------------------------------------------------------------- endif else begin xrange = [-1, 1] yrange = [ 0, 1000] ierror = 1B endelse ; PLOT THE ***TRANSFER FUNCTIONS*** AS A FUNCTION OF WAVELENGTHS ;------------------------------------------------------------------------------- loadct, 0, /silent plot, [5.0, 14.0], yrange, xstyle=5, ystyle=4, $ positi=[13800, 4350, 17800, 7350], /device, /noerase, /nodata for i=0,n_elements(wavele)-1 do polyfill, wavele[i]+[-1,+1,+1,-1,-1]*wavdel, $ [ 1, 1, 0, 0, 1]*!y.crange[0]+[ 0, 0, 1, 1, 0]*!y.crange[1], color=220 if not ierror then begin loadct, 33, /silent ; LOOP OVER THE INDIVIDUAL OBJECTS TO PLOT THEM ;----------------------------------------------------------------------- for i=0,n_elements(a2jdat)-1 do begin tmpvar = 165 + (1.0*i/(n_elements(a2jdat)-1))^1.2 * (254-164) oplot, a2jdat[i].wavele[n_elements(wavele):*], $ a2jdat[i].(tagidx)[n_elements(wavele):*], color=tmpvar if a2jdat[i].qualit eq 2 then badstr = ['', '*'] $ else badstr = ['', ''] if flgdat[i] eq 2 then badstr += ['(', ')'] $ else badstr += [' ', ''] xyouts, 5.4, (0.93-i/20.)*!y.crange[1], badstr[0]+ $ strtrim(a2jdat[i].object,2)+badstr[1], color=tmpvar, $ charsize=0.7*!p.charsize endfor ; PLOT THE AVERAGE AND THE ERROR OF THE RESPECTIVE CONVERSION FACTOR ;----------------------------------------------------------------------- loadct, 39, /silent oplot, a2jave.wavele[n_elements(wavele):*], $ a2jave.(tagidx )[n_elements(wavele):*], thick=2, linestyle=5 oplot, a2jave.wavele[n_elements(wavele):*], $ a2jave.(tagidx+1)[n_elements(wavele):*], thick=2, linestyle=1 endif else begin xyouts, 15800, 2000, 'No calib data found.', $ charsize=1.8*!p.charsize, align=0.5, /device endelse if !p.font gt -1 then strvar = ['!9l!X [!9m!Xm]'] $ else strvar = ['!4k!X [!4l!Xm]'] plot, [5.0, 14.0], yrange, xstyle=1, positi=[13800, 4350, 17800, 7350], $ xtitle='wavelength '+strvar, /device, /noerase, /nodata xyouts, 13100, 5850, 'ADU / Jy', align=0.5, orientation=90, /device xyouts, 15800, 7450, 'conversion factors', align=0.5, $ charsize=1.5*!p.charsize, /device ; PLOT THE ***TRANSFER FUNCTIONS*** AS A FUNCTION OF TIME ;------------------------------------------------------------------------------- plot, xrange, yrange, xstyle=4, ystyle=4, positi=[13800, 550, 17800, 3550], $ /device, /noerase, /nodata if not ierror then begin ; DECOMPOSE THE QUALITY FLAGS INTO AN ARRAY ;----------------------------------------------------------------------- qualit = reverse(byte((string(a2jdat.qualit, format='(B08)')))-48B) ; PLOT THE CALIBRATOR NAMES ;----------------------------------------------------------------------- loadct, 0, /silent for i=0,n_elements(a2jdat)-1 do begin if qualit[5,i] then badstr = ['', '*'] $ else badstr = ['', ''] if flgdat[i] eq 2 then badstr += ['(', ')'] $ else badstr += [' ', ''] xyouts, a2jdat[i].mjdobs+0.01*(!x.crange[1]-!x.crange[0]), $ 0.05*!y.crange[1], badstr[0]+ $ strtrim(a2jdat[i].object,2)+badstr[1], color=210, $ orientat=90, charsize=0.5*!p.charsize endfor ; LOOP OVER THE WAVELENGTHS ;----------------------------------------------------------------------- loadct, 39, /silent for i=0,n_elements(wavele)-1 do begin ; PLOT ALL THE GOOD CALIBRATORS WITH FILLED CIRCLES ;----------------------------------------------------------------------- select = where(flgdat eq 0) if select[0] ge 0 then begin usersym, cos(findgen(17) * (!pi*2/16.)), $ sin(findgen(17) * (!pi*2/16.)), /fill plots, a2jdat[select].mjdobs, a2jdat[select].corflx[i], $ color=colour[i], psym=8, symsize=0.5 errplot, a2jdat[select].mjdobs, color=colour[i], width=0, $ a2jdat[select].corflx[i]-a2jdat[select].corflxerr[i], $ a2jdat[select].corflx[i]+a2jdat[select].corflxerr[i] endif ; PLOT ALL THE REJECTED CALIBRATORS WITH OPEN BOXES ;----------------------------------------------------------------------- select = where(flgdat eq 2) if select[0] ge 0 then begin usersym, cos(findgen(17) * (!pi*2/16.)), $ sin(findgen(17) * (!pi*2/16.)) plots, a2jdat[select].mjdobs, a2jdat[select].corflx[i], $ color=colour[i], psym=8, symsize=0.5 errplot, a2jdat[select].mjdobs, color=colour[i], width=0, $ a2jdat[select].corflx[i]-a2jdat[select].corflxerr[i], $ a2jdat[select].corflx[i]+a2jdat[select].corflxerr[i] endif ; PLOT THE AVERAGE TRANSFER FUNCTION AND THE LOCATION OF CURRENT DATA ;----------------------------------------------------------------------- oplot, !x.crange, a2jave[0].corflx[i]*[1,1], color=colour[i], $ linestyle=1 if calidx ge 0 then tmpvar = a2jdat[calidx].corflx[i] $ else tmpvar = a2jave[0].corflx[i] plots, 0, tmpvar, psym=7, symsize=0.5, thick=2, color=colour[i] errplot, 0, tmpvar-a2jave[0].corflxerr[i], $ tmpvar+a2jave[0].corflxerr[i], $ color=colour[i], thick=2 ; MARK THE CALIBRATOR USED FOR THE CURRENT DATA ;----------------------------------------------------------------------- if calidx ge 0 then begin plots, a2jdat[calidx].mjdobs, a2jdat[calidx].corflx[i], $ color=colour[i], psym=6, symsize=0.6 endif ; END OF LOOP OVER THE WAVELENGTHS ;----------------------------------------------------------------------- endfor endif else begin xyouts, 15800, 2000, 'No calib data found.', $ charsize=1.8*!p.charsize, align=0.5, /device endelse if !p.font gt -1 then strvar = ['!9D!X'] $ else strvar = ['!4D!X'] plot, xrange, yrange, positi=[13800, 550, 17800, 3550], $ xtitle='time '+strvar+'t [h]', /device, /noerase, /nodata xyouts, 13100, 2050, 'ADU / Jy', align=0.5, orientation=90, /device ; TRY TO LOAD A FILE WITH THE UV COORDINATES ;------------------------------------------------------------------------------- tmpvar = where(file_test(ewstag+filnam[1:2,*]) eq 1) if tmpvar[0] ge 0 then begin tmpvar = ewstag+(filnam[1:2,*])[tmpvar[0]] tmpdat = oirgetvis(tmpvar, ierr=ierror) endif else ierror = 1B ; CHECK IF THE DATA WAS LOADED SUCCESSFULLY ;------------------------------------------------------------------------------- if ~ ierror then begin ; THEN GET THE UV COORDINATES AND DETRMINE THE SCALING OF THE PLOT ;----------------------------------------------------------------------- if (tmpdat.ucoord lt 0.0) then begin uvcoor = [-tmpdat.ucoord, -tmpdat.vcoord] endif else begin uvcoor = [+tmpdat.ucoord, +tmpdat.vcoord] endelse uvsize = 1.05*max(sqrt(total(uvcoor^2))) ; THEN GET THE NAMES OF THE STATIONS ;----------------------------------------------------------------------- strvar = midi_getkey(tmpvar, ['ISS CONF STATION1', 'ISS CONF STATION2']) strvar = strtrim(strvar.(1),2) + '-' + strtrim(strvar.(2),2) endif else uvsize = 1.0 ; PLOT A SMALL UV PLANE WITH THE LOCATION OF THE MEASUREMENTS ;------------------------------------------------------------------------------- usersym, cos(findgen(17) * (!pi*2/16.)), sin(findgen(17) * (!pi*2/16.)), /fill !x.ticklen = 0.026 !y.ticklen = 0.052 loadct, 0, /silent plot, [uvsize, 0.0001], [-uvsize, uvsize], xstyle=5, ystyle=5, $ xrange=[uvsize, 0.0001], xmargin=[0,0], ymargin=[0,0], $ position=[9400, 10000, 12000, 15200], /device, /nodata, /noerase oplot, -sin(findgen(101)/50. * !pi) * fix(uvsize/10-0.5) * 5, $ +cos(findgen(101)/50. * !pi) * fix(uvsize/10-0.5) * 5, color=175 oplot, -sin(findgen(101)/50. * !pi) * fix(uvsize/10-0.5) * 10, $ +cos(findgen(101)/50. * !pi) * fix(uvsize/10-0.5) * 10, color=175 loadct, 39, /silent if uvsize ne 1.0 then begin oplot, [0,uvcoor[0]], [0,uvcoor[1]], color=80 plots, uvcoor[0], uvcoor[1], psym=8, color=80, symsize=1.0 xyouts, uvcoor[0], uvcoor[1]+(0.12*uvcoor[1]/abs(uvcoor[1]) $ -0.03)*uvsize, strvar, color=80, align=0.5 endif axis, 0 , 0, /xaxis, xstyle=1, xtickinterval=fix(uvsize/10-0.5)*5, $ xtickname=replicate(' ', 10), ticklen=0.03 axis, 0 , 0, /yaxis, ystyle=1, ytickinterval=fix(uvsize/10-0.5)*5, $ ticklen=0.06 xyouts, 0.07 * uvsize, 1.03 * uvsize, 'V' xyouts, 1.07 * uvsize, -0.07 * uvsize, 'U' xyouts, 10700, 15650, 'UV plane', charsize=1.5*!p.charsize, $ align=0.5, orientation=0, /device ; CLOSE THE PLOTTING DEVICE AND RESET THE PLOT PARAMETERS ;------------------------------------------------------------------------------- !x.ticklen = 0.0 !y.ticklen = 0.0 device, /close set_plot, pltset.d_old !p = pltset.p END ;+ ; NAME: ; MIDI_PLTRAD ; ; PURPOSE: ; Plot the correlated flux or visibility as a function of baseline length. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; midi_pltrad, caldat [, wavele] [, wavdel] [, /corflx] [, /visamp] ; [, xrange=xrange] [, yrange=yrange] [, /ylog] ; [, titles=titles] [, colour=colour] ; [, psfile=psfile] ; ; REQUIRED INPUTS: ; caldat MIDI result structure ; ; OPTIONAL INPUTS: ; wavele wavelengths for which the data is to be plotted ; wavdel half width of the wavelength range around WAVELE ; ; KEYWORDS: ; corflx make a plot using the correlated fluxes (default) ; visamp make a plot using the visibilities ; xrange range of the x axis (baseline length in m) ; yrange range of the y axis (correlated flux or visibility) ; ylog use logarithmic scaling for the y-axis ; titles string containing a custom title for the plot ; colour colour numbers for the different wavelengths ; psfile selects output to PS and specifies the file name ; ; OUTPUTS: ; plot of the correlated flux of visibility vs. baseline length ; ; DESCRIPTION AND EXAMPLE: ; This programme creates a plot of either the correlated flux or the ; visibility as a function of the projected baseline length. By default ; a plot for the correlated flux is created; visibilities are plotted ; with the keyword VISAMP set. The programme uses a MIDI result struc- ; ture as an input: CALDAT. For every element in CALDAT, the correlated ; flux or visibility is averaged around the wavelength specified in ; WAVELE (default is 12 micron) in a range with the half width given by ; WAVDEL (default value is 0.2 micron). The same is done for the errors ; of the correlated fluxes or visibilities. These averages are then ; plotted as a function of the projected baseline lengths, which are ; calculated from the uv-coordinates given in CALDAT. If WAVELE is an ; array with several wavelengths, then the values are calculated and ; plotted for multiple wavelengths. Specific colours for the different ; wavelengths can be specified with the parameter COLOUR. In the ; simplest case, the programme is called by: ; midi_pltrad, caldat ; ; CALLED BY: ; none ; ; CALLING: ; setplot ; ; MODIFICATION HISTORY: ; 2010-10-05 Written by Konrad R. Tristram ; PRO MIDI_PLTRAD, CALDAT, WAVELE, WAVDEL, CORFLX=CORFLX, VISAMP=VISAMP, $ XRANGE=XRANGE, YRANGE=YRANGE, YLOG=YLOG, TITLES=TITLES, $ ERRWID=ERRWID, COLOUR=COLOUR, PSFILE=PSFILE ; SET DEFAULT VALUE FOR THE WAVELENGTH RANGE TO BE PLOTTED ;------------------------------------------------------------------------------- if n_elements(wavele) lt 1 then wavele = 12.0 numwav = n_elements(wavele) if n_elements(wavdel) ne 1 then wavdel = 0.2 ; SELECT IF PLOTTING CORRELATED FLUX OR VISIBILITY ;------------------------------------------------------------------------------- if keyword_set(visamp) then begin tagidx = where(tag_names(caldat) eq 'VISAMP') ytitle = 'visibility V' legend = 'visibility at ' endif else begin tagidx = where(tag_names(caldat) eq 'CORFLX') ytitle = 'correlated flux F!Dcor!N [Jy]' legend = 'corr. flux at ' endelse ; CALCULATE AN AUTOMATIC COLOUR SCALE ;------------------------------------------------------------------------------- if n_elements(colour) ne numwav then begin colour = 50 + 200*findgen(numwav)/((numwav-1) > 1) endif ; CREATE THE VARIABLES HOLDING THE VALUES TO PLOTTED ;------------------------------------------------------------------------------- x_vals = fltarr(numwav, n_elements(caldat)) y_vals = x_vals y_errs = x_vals ; SAVE THE VALUES TO BE PLOTTED IN THE VARIABLES (NEEDED FOR AUTOMATIC SCALING) ;------------------------------------------------------------------------------- for i=0,numwav-1 do begin for j=0L,n_elements(caldat)-1 do begin select = where((caldat[j].wavele ge (wavele[i]-wavdel)) and $ (caldat[j].wavele le (wavele[i]+wavdel))) if select[0] ge 0 then begin x_vals[i,j] = sqrt(caldat[j].ucoord^2 + $ caldat[j].vcoord^2) y_vals[i,j] = mean(caldat[j].(tagidx )[select]) y_errs[i,j] = mean(caldat[j].(tagidx+1)[select]) endif else begin x_vals[i,j] = !values.f_nan y_vals[i,j] = !values.f_nan y_errs[i,j] = !values.f_nan endelse endfor endfor ; CALCULATE THE SCALING OF THE PLOT ;------------------------------------------------------------------------------- if n_elements(xrange) ne 2 then xrange=[0.0, max(x_vals)] else xstyle=1 if n_elements(yrange) ne 2 then yrange=[1E-3, max(y_vals+y_errs)] else ystyle=1 ; CALCULATE THE SCALING OF THE PLOT ;------------------------------------------------------------------------------- if keyword_set(psfile) then symsiz=0.75 else symsiz=1.25 ; OPEN THE PLOT AND PLOT THE AXIS ;------------------------------------------------------------------------------- loadct, 0, /silent setplot, pltset, pltpos, psfile=psfile plot, xrange, yrange, xtitle='projected baseline length BL [m]', $ ytitle=ytitle, xstyle=xstyle, ystyle=ystyle, position=pltpos, $ title=titles, ylog=ylog, /device, /nodata, /noerase ; OVERPLOT VISIBILITY=1 FOR THE VISIBILITY PLOT ;------------------------------------------------------------------------------- if keyword_set(visamp) then oplot, [-1D5, +1D5], [1,1], linestyle=2, color=200 ; PLOT THE DATA AND LEGEND ;------------------------------------------------------------------------------- loadct, 39, /silent if !p.font gt -1 then strvar = ['" !9m!Xm"'] $ else strvar = ['" !4l!Xm"'] for i=0,numwav-1 do begin oplot, x_vals[i,*], y_vals[i,*], psym=8, symsize=symsiz, color=colour[i] errplot, x_vals[i,*], y_vals[i,*]-y_errs[i,*], width=errwid, $ y_vals[i,*]+y_errs[i,*], color=colour[i] plots, pltpos[0]+ 0.65*(pltpos[2]-pltpos[0]), $ pltpos[1]+(0.94-i*0.05)*(pltpos[3]-pltpos[1]), $ psym=8, color=colour[i], symsize=symsiz, /device xyouts, pltpos[0]+ 0.67*(pltpos[2]-pltpos[0]), $ pltpos[1]+(0.93-i*0.05)*(pltpos[3]-pltpos[1]), $ legend+string(wavele[i], format='(F4.1, '+strvar+')'), /device endfor ; CLOSE THE PLOT AND RESET PLOT PARAMETERS ;------------------------------------------------------------------------------- setplot, pltset, /close END ;+ ; NAME: ; MIDI_PLTUVP ; ; PURPOSE: ; Plot points in a UV plane with values from a MIDI result struct. ; ; CATEGORY: ; MIDI data tools. ; ; CALLING SEQUENCE: ; midi_pltuvp, caldat [, wavele] [, tagnam=tagnam] [, zrange=zrange] ; [, declin=declin] [, blname=blname] [, minalt=minalt] ; [, pblran=pblran] [, titles=titles] [, sqroot=sqroot] ; [, barlab=barlab] [, /noline] [, labels=labels] ; [, /oldwin] [, psfile=psfile] ; ; REQUIRED INPUTS: ; caldat MIDI result structure ; ; OPTIONAL INPUTS: ; wavele wavelength for which the data is to be plotted ; ; KEYWORDS: ; tagnam name of the tag from the struct to be plotted ; zrange range of the colour scaling ; declin declination of the target ; blname baseline configurations (as string array or indices) ; minalt minimum altitude the object must have ; pblran range of baseline lengths ; titles title string for the plot ; sqroot use square root scaling of the color map ; barlab label of the colour bar ; noline don't plot the cross hairs ; labels label the UV points with numbers and/or values ; oldwin reuse the old plotting window (only for screen output) ; pltsiz size of the plot (A&A single-column, 10cm, etc.) ; psfile name of the PS file for PS output ; ; OUTPUTS: ; a PS file or plot on the screen ; ; DESCRIPTION AND EXAMPLE: ; This programme makes a plot of the data in a MIDI result struct on ; the UV plane. For every UV point, a point, colour coded corresponding ; to its value, is plotted in the UV plane. The programme first uses ; VLTI_PLTUVP to create the basic plot of the UV plane and it offers ; many of the options of that routine for customising the plot. Then ; the data values are overplotted. The wavelength for which the values ; are to be plotted can be specified using WAVELE: the first element in ; WAVELE specifies the central wavelength and the second element the ; half width of the wavelength band to be averaged. The default is ; WAVELE = [12.0, 0.2], that is all bins with 11.8 < wavele < 12.2 are ; averaged. By default the correlated flux (i.e. CALDAT.CORFLX) is ; plotted; a different quantity can be chosen by specifying the corres- ; ponding tag name in the keyword TAGNAM. By default the scaling of the ; values is linear. Optionally, a square root colour scale can be used ; by setting the keyword SQROOT. ; The simplest call of the function is by: ; midi_pltuvp, caldat ; To plot the visibility from 11 to 12 micron into a postscript file ; for an object at DEC = -65 and the 'E0-G0' and 'U2-U4' baselines as ; well as with customised labels and titles the call is: ; midi_pltuvp, mydata, tagnam='visamp', declin=-65, minalt=25., $ ; blname=['E0-G0', 'U2-U4'], pblran=[-110,110], $ ; titles='My data', barlab='visibility V', psfile='myplot' ; ; CALLED BY: ; none ; ; CALLING: ; vlti_pltuvp ; plotcbar ; ; MODIFICATION HISTORY: ; 2010-12-03 Written by Konrad R. Tristram ; 2012-03-29 K. Tristram: Added keyword PLTSIZ to set size of the plot. ; PRO MIDI_PLTUVP, CALDAT, WAVELE, TAGNAM=TAGNAM, ZRANGE=ZRANGE, DECLIN=DECLIN, $ BLNAME=BLNAME, MINALT=MINALT, PBLRAN=PBLRAN, TITLES=TITLES, $ SQROOT=SQROOT, BARLAB=BARLAB, NOLINE=NOLINE, LABELS=LABELS, $ OLDWIN=OLDWIN, PLTSIZ=PLTSIZ, PSFILE=PSFILE ; GET THE NUMBER OF UV POINTS ;------------------------------------------------------------------------------- numele = n_elements(caldat) ; CHECK OR SET THE WAVELENGTHS AND DETERMINE THEIR INDICES ;------------------------------------------------------------------------------- case n_elements(wavele) of 1: wavele = [float(wavele) > 7.5 < 13.5, 0.2] 2: wavele = float(wavele) else: wavele = [12.0, 0.2] endcase welmnt = where((caldat[0].wavele le wavele[0]+wavele[1]) and $ (caldat[0].wavele ge wavele[0]-wavele[1])) ; CHECK IF A CORRECT TAG NAME IS GIVEN ;------------------------------------------------------------------------------- if ~ keyword_set(tagnam) then tagnam = 'corflx' tagidx = where(tagnam eq strlowcase(tag_names(caldat))) if tagidx[0] lt 0 then tagidx = where('corflx' eq strlowcase(tag_names(caldat))) ; GET THE DATA TO BE PLOTTED FROM THE RESULT STRUCT ;------------------------------------------------------------------------------- z_data = reform(caldat.(tagidx)[welmnt], n_elements(welmnt), numele) z_data = total(z_data, 1)/n_elements(welmnt) ; DETRMINE ZRANGE AND SCALE THE DATA FOR PLOTTING ;------------------------------------------------------------------------------- if n_elements(zrange) ne 2 then zrange = [min(z_data), max(z_data)] z_colo = findgen(256) if keyword_set(sqroot) then begin z_plot = sqrt((z_data-zrange[0])/(zrange[1]-zrange[0]))*254 > 0 < 254 z_valu = (z_colo/255)^2 * (zrange[1]-zrange[0]) + zrange[0] z_colo = sqrt(z_colo/255.)*255. endif else begin z_plot = ((z_data-zrange[0])/(zrange[1]-zrange[0]))*254 > 0 < 254 z_valu = (z_colo/255) * (zrange[1]-zrange[0]) + zrange[0] endelse ; SET THE PLOTTING RANGE ;------------------------------------------------------------------------------- if n_elements(pblran) ne 2 then begin pblran = max(abs([caldat.ucoord, caldat.vcoord])) tmpvar = 10.^(round(alog10(pblran) + 0.4) - 1.0) pblran = [-1.,+1.] * round(1.4 * pblran / tmpvar) * tmpvar endif ; OPEN THE UV PLANE PLOT ;------------------------------------------------------------------------------- vlti_pltuvp, pltset, declin, blname=blname, minalt=minalt, pblran=pblran, $ titles=titles, barpos=barpos, noline=noline, nolabel=nolabel, $ oldwin=oldwin, pltsiz=pltsiz, psfile=psfile ; PLOT THE LABELS FOR THE UV POINTS ;------------------------------------------------------------------------------- loadct, 0, /silent if keyword_set(psfile) then chrsiz = 0.25 else chrsiz = 0.65 if keyword_set(labels) then for i=0L,numele-1 do begin case labels of 1: tmpvar = string(z_data[i], format='(F5.2)') 2: tmpvar = string(i, format='(I3)') else: tmpvar = string(i, z_data[i], format='(I3, ":", F5.2)') endcase xyouts, -caldat[i].ucoord-max(pblran)/50., -caldat[i].vcoord, tmpvar, $ charsize=chrsiz, color=200, orientation=60 endfor ; PLOT THE UV POINTS ;------------------------------------------------------------------------------- loadct, 39, /silent if keyword_set(psfile) then chrsiz = 0.75 else chrsiz = 1.0 for i=0L,numele-1 do begin oplot, caldat[i].ucoord*[-1.,+1.], caldat[i].vcoord*[-1.,+1.], $ psym=8, color=z_plot[i], symsize=chrsiz endfor ; PLOT THE COLOUR BAR ;------------------------------------------------------------------------------- if ~ keyword_set(barlab) then barlab = tagnam plotcbar, z_valu, z_colo, barpos, barlab=barlab ; majval=[0.5, 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11.], $ ; minval=[findgen(10)/10., findgen(20)/5.+1.0, findgen(15)/2.+5.0], $ ; format='(F4.1)', /endlab ; CLOSE THE UV PLANE PLOT ;------------------------------------------------------------------------------- vlti_pltuvp, pltset, /close END