; MIDI_PREPHT - Prefix the photometry to a MIDI result structure. ; outdat = midi_prepht(orgdat, flxsel, /noweight) ; ; MIDI_EXTVIS - Extract visibilities from a model at specified UV coordinates. ; result = midi_extvis(uvplan, uvcoor, rotang, /uvplot) ; ; MIDI_DATCOM - Interpolate MIDI data at certain wavelengths. ; result = midi_datcom(caldat, wavele, smtval, wavwid=0.0) ; ; MIDI_DPHASE - Determine the differential phase from input phase data. ; result = midi_dphase(mdlvis, method, wavran=[8,13], smtval=3) ; ; MIDI_PLTCMP - Plot a comparison of measured data and modelled data. ; midi_pltcmp, datvis, mdlvis, mdlvi1, wavran=[8,13], /normal ; ;+ ; NAME: ; MIDI_PREPHT ; ; PURPOSE: ; Prefix the photometry to a MIDI result structure. ; ; CATEGORY: ; MIDI model comparison. ; ; CALLING SEQUENCE: ; outdat = midi_prepht(orgdat [, flxsel] [, /noweight]) ; ; REQUIRED INPUTS: ; orgdat input array containing the MIDI results ; ; OPTIONAL INPUTS: ; flxsel array containing the indices of elements to consider ; ; KEYWORDS: ; noweight use weighted averaging instead of unweighted averaging ; totflx use total flux instead of masked flux as correlated flux ; ; OUTPUTS: ; outdat output array in the same format as ORGDAT ; ; DESCRIPTION AND EXAMPLE: ; This function adds a new array element at the beginning of the input ; structure and defines it as a baseline with a baseline length BL = 0, ; The input data must be an array of structure MIDI result, that is ; the result of MIDI measurements (e.g. either "MIDI_PRISM" or ; "MIDI_GRISM"). The arithmetic mean of the selected total fluxes is ; saved as the correlated flux of the new element. The mean can be ; either calculated by weighting the individual measurements with the ; inverse square of their errors or by using no weighting of the ; measurements (if the keyword NOWEIGHT is set). In the first case, ; the error of the averaged total flux is determined from the errors ; of the individual measurements; in the second case the error of the ; averaged total flux is determined by the standard deviation of the ; individual measurements (which of course only makes sense for a ; certain number of measurements). The elements to be used for the ; determination of the averaged total flux can be specified in FLXSEL. ; If this variable is not given, all elements are used. Example: ; result = prefphot(mydata) ; ; CALLED BY: ; midi_2bbpre ; midi_3bbpre (and other fit preparing routines) ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2007-02-05 Adapted from earlier test_5.pro by Konrad R. Tristram ; 2008-09-17 Konrad R. Tristram: Convert to the new MIDI result format. ; 2009-07-20 K. Tristram: Renamed function and added keyword NOWEIGHT. ; 2012-11-12 K. Tristram: Masked -> total flux, add keyword TOTFLX. ; FUNCTION MIDI_PREPHT, ORGDAT, FLXSEL, NOWEIGHT=NOWEIGHT, TOTFLX=TOTFLX ; CHECK IF FLXSEL IS PRESENT OTHERWISE SELECT ALL ELEMENTS ; ----------------------------------------------------------------------------- if (n_elements(flxsel) le 0) then flxsel = indgen(n_elements(orgdat)) flxsel = long(flxsel) > 0 < (n_elements(orgdat)-1) ; GET THE NUMBER OF WAVELENGTHS AND NUMBER OF ELEMENTS IN THE STRUCTURE ; ----------------------------------------------------------------------------- numwav = n_elements(orgdat[0].wavele) numele = n_elements(orgdat) ; CREATE THE STRUCTURE WHICH WILL HOLD THE RESULT ; ----------------------------------------------------------------------------- resdat = [orgdat[0], reform(orgdat, numele)] ; SET STRAIGHT FORWARD VALUES FOR THE ZERO BASELINE ; ----------------------------------------------------------------------------- resdat[0].wavele = orgdat[0].wavele resdat[0].baseli = 'total' resdat[0].ucoord = 0.0 resdat[0].vcoord = 0.0 resdat[0].pbline = 0.0 resdat[0].pangle = 0.0 resdat[0].visamp = 1.0 resdat[0].visamperr = 0.001 resdat[0].visphi = 0.0 resdat[0].visphierr = 0.001 ; CHECK IF ERRORS ARE SUPPOSED TO BE USED AS WEIGHTS ; ----------------------------------------------------------------------------- if keyword_set(noweight) then begin ; USE THE STANDARD DEVIATION OF THE TOTAL AND MASKED FLUXES AS ERRORS ;---------------------------------------------------------------------- for i=0,numwav-1 do begin resdat[0].totflxerr[i] = stddev(orgdat[flxsel].totflx[i]) / $ sqrt(n_elements(flxsel)) resdat[0].mskflxerr[i] = stddev(orgdat[flxsel].mskflx[i]) / $ sqrt(n_elements(flxsel)) endfor ; USE UNIFORM WEIGHTING ;---------------------------------------------------------------------- totwei = make_array(numwav, n_elements(flxsel), value=1.0) mskwei = make_array(numwav, n_elements(flxsel), value=1.0) endif else begin ; CALCULATE THE WEIGHTED ERRORS OF THE AVERAGED FLUXES ;---------------------------------------------------------------------- resdat[0].totflxerr = sqrt(1./total(1./orgdat[flxsel].totflxerr^2, 2)) resdat[0].mskflxerr = sqrt(1./total(1./orgdat[flxsel].mskflxerr^2, 2)) ; USE THE INVERSE SQUARES OF THE ERRORS AS WEIGHTS FOR THE AVERAGING ;---------------------------------------------------------------------- totwei = 1.0 / orgdat[flxsel].totflxerr^2 mskwei = 1.0 / orgdat[flxsel].totflxerr^2 endelse ; CALCULATE THE (WEIGHTED) ARITHMETIC MEAN OF THE TOTAL AND MASKED FLUXES ; ----------------------------------------------------------------------------- resdat[0].totflx = total(orgdat[flxsel].totflx * totwei, 2) / total(totwei, 2) resdat[0].mskflx = total(orgdat[flxsel].mskflx * mskwei, 2) / total(mskwei, 2) ; COPY THE FLUX VALUES AND ERRORS OF THE TOTAL FLUX TO THE CORRELATED FLUX ; ----------------------------------------------------------------------------- if keyword_set(totflx) then begin resdat[0].corflx = resdat[0].totflx resdat[0].corflxerr = resdat[0].totflxerr endif else begin resdat[0].corflx = resdat[0].mskflx resdat[0].corflxerr = resdat[0].mskflxerr endelse ; RETURN THE RESULT ; ----------------------------------------------------------------------------- return, resdat END ;+ ; NAME: ; MIDI_EXTVIS ; ; PURPOSE: ; Extract visibilities from a model at specified UV coordinates. ; ; CATEGORY: ; MIDI model comparison. ; ; CALLING SEQUENCE: ; result = midi_extvis(uvplan, uvcoor [, rotang] [, /uvplot]) ; ; REQUIRED INPUTS: ; uvplan structure containing the fourier transform of the model ; uvcoor structure containing the UV coordinates to be extracted ; ; OPTIONAL INPUTS: ; rotang additional model rotation angle, counted anticlockwise ; ; KEYWORDS: ; uvplot make plot of the uv plane with location of UV coordinates ; direct do the interpolation directly in the complex UV plane ; ; OUTPUTS: ; result structure containing the extracted visibilities ; ; DESCRIPTION AND EXAMPLE: ; This function extracts the visbilities at certain UV points from the ; UV plane of a model. The input model, UVPLAN, must be a UV plane of ; the model, as generated for example by MODL_M2VISI. UVPLAN also ; determines the wavelengths for which the visibilities are determined. ; The UV points are selected by the UV coordinates given in UVCOOR. ; UVCOOR can either be measured data (e.g. of type MIDI_PRISM) or a set ; of UV tracks, e.g. produced by 'uvcoor = vlti_blproj(-65.0, ['U2-U3', ; 'U3-U4'], hrange=[-4.0,+4.0])'. The programme then calculates the x ; and y coordinates in the array of the Fourier transform for the given ; wavelength, pixel scale and baseline information. Then the correlated ; flux and the visibility are interpolated at these coordinates. By ; default, a subarray is cut out from the complex UV plane at every UV ; point and the visibility and phase are calculated from it. The inter- ; polation is then carried out in this subarray. A simple algorithm ; ensures that there are no errors in the interpolation of the phase ; caused by a wrapping of the phase. If the keyword DIRECT is set, then ; the complex visibility is directly interpolated in the complex UV ; plain and only then are the visibility and phase calculated from it. ; This second method is, however, not very accurate for strong gradients ; in the phase. The total flux is determined by the maximum of the ; Fourier transform. The resulting structure is similar to the one of ; measured data, except that the number of wavelengths in every UV point ; is different and that all errors are set to 0.0. Optionally a plot of ; the UV plane with the baseline positions are generated. ; The function is invoked by: result = midi_extvis(uvplan, uvcoora) ; ; CALLED BY: ; midi_2bbdat ; midi_2bbcmp ; midi_3bbdat ; midi_3bbcmp ; ; CALLING: ; modl_m2plot ; ; MODIFICATION HISTORY: ; 2006-05-03 Written by Konrad R. W. Tristram ; 2006-06-26 K. Tristram: Corrected direction of the RA axis in the plot. ; 2006-07-03 K. Tristram: Also transfer u and v coordinates to result. ; 2007-02-07 K. Tristram: Added the extraction of phases. ; 2008-09-17 K. Tristram: Adapted to use the new MIDI result structure. ; 2009-07-13 K. Tristram: Renamed function and various other changes. ; 2009-08-17 K. Tristram: Correct random phases for very low visibility. ; 2010-04-03 K. Tristram: Correct subscript error for baseline copying. ; 2010-09-23 K. Tristram: Make result compatible to standard MIDI struct. ; FUNCTION MIDI_EXTVIS, UVPLAN, UVCOOR, ROTANG, UVPLOT=UVPLOT, DIRECT=DIRECT ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- if (n_elements(rotang) eq 0) then rotang = 0.0 else rotang = float(rotang) uvpelm = [n_elements(uvcoor[0].pbline), n_elements(uvcoor)] wavelm = n_elements(uvplan) overpl = 0 ; CREATE THE STRUCTURE HOLDING THE RESULT ;------------------------------------------------------------------------------- result = {identi:'', object:'', pbline:0.0, pangle:0.0, $ ucoord:0.0, vcoord:0.0, baseli:'', mjdobs:0.0D, $ wavele:uvplan.wavele, $ mskflx:fltarr(wavelm), mskflxerr:fltarr(wavelm), $ totflx:fltarr(wavelm), totflxerr:fltarr(wavelm), $ corflx:fltarr(wavelm), corflxerr:fltarr(wavelm), $ visamp:fltarr(wavelm), visamperr:fltarr(wavelm), $ visphi:fltarr(wavelm), visphierr:fltarr(wavelm), $ commen:''} result = replicate(result, uvpelm) ; COPY BASELINE INFORMATION ;------------------------------------------------------------------------------- i = n_elements(uvpelm)-1 while (uvpelm[i] eq 1) and (i gt 0) do i -= 1 uvpelm = uvpelm[0:i] result.pbline = reform(uvcoor.pbline, uvpelm) result.pangle = reform(uvcoor.pangle, uvpelm) result.ucoord = reform(uvcoor.ucoord, uvpelm) result.vcoord = reform(uvcoor.vcoord, uvpelm) for i=0,n_elements(uvcoor)-1 do result[*,i].baseli = uvcoor[i].baseli ; GET OBJECT DESCRIPTION IF IT IS DEFINED IN THE MODEL STRUCTURE ;------------------------------------------------------------------------------- if total(tag_names(uvplan) eq 'OBJECT') then begin object = uvplan.object ; IF THERE ARE DIFFERENT DESCRIPTIONS THEN CONCATINATE THEM ;----------------------------------------------------------------------- if ~ array_equal(object, shift(object, 1)) then begin result.object = strjoin(object, '; ') ; OTHERWISE SAVE THE FIRST OBJECT DESCRIPTION INTO THE RESULT ;----------------------------------------------------------------------- endif else result.object = object[0] endif ; LOOP OVER THE ELEMENTS IN THE MODEL ;------------------------------------------------------------------------------- for i=0,wavelm-1 do begin ; PREPARE THE CALCULATION OF THE POSITIONS IN THE UV-PLANE ;----------------------------------------------------------------------- arrsiz = size(uvplan[0].fftdat) factor = !pi / (uvplan[i].wavele * 1.0e-6) / (3600. * 180. * 1000.) * $ arrsiz[1] * uvplan[i].pscale ; USE THE BASELINE IN POLAR COORDINATES FOR THE CALCULATION ;----------------------------------------------------------------------- radius = factor * result.pbline x = arrsiz[1]/2. - sin((result.pangle - rotang)/180.*!pi) * radius y = arrsiz[2]/2. + cos((result.pangle - rotang)/180.*!pi) * radius ; USE THE U AND V COORDINATES FOR THE CALCULATION ;----------------------------------------------------------------------- ;x = arrsiz[1]/2. - factor * result.ucoord ;y = arrsiz[2]/2. + factor * result.vcoord ; DETERMINE THE TOTAL FLUX FROM THE MAXIMUM OF THE FOURIER TRANSFORM ;----------------------------------------------------------------------- totflx = abs(max(uvplan[i].fftdat)) result.mskflx[i] = totflx result.totflx[i] = totflx ; IF DIRECT IS SET, THEN DIRECTLY INTERPOLATE IN THE COMPLEX VALUES ;----------------------------------------------------------------------- if keyword_set(direct) then begin ; INTERPOLATE THE COMPLEX VISIBILITY AT THE UV COORDINATES ;--------------------------------------------------------------- comvis = interpolate(uvplan[i].fftdat, x, y, cubic=-0.5) ; CALCULATE THE VISIBILITY AND PHASE FROM THE COMPLEX VALUE ;--------------------------------------------------------------- result.corflx[i] = abs(comvis) result.visphi[i] = atan(comvis, /phase) * 180. / !pi ; OTHERWISE CUT OUT SUBARRAY AND CALCULATE THE REAL NUMBERS FROM IT ;----------------------------------------------------------------------- endif else for j=0,n_elements(x)-1 do begin ; SET SIZE OF THE SUBARRAY AND CALCULATE THE NEW COORDINATES ;--------------------------------------------------------------- ss = 2 xx = ss + (x[j] mod 1) yy = ss + (y[j] mod 1) ; GET THE SUBARRAY AND INTERPOLATE THE CORRELATED FLUX ;--------------------------------------------------------------- subarr = uvplan[i].fftdat[x[j]-ss:x[j]+ss, y[j]-ss:y[j]+ss] result[j].corflx[i] = interpolate(abs(subarr), xx, yy, cub=-0.5) ; CALCULATE THE PHASES FROM THE COMPLEX VISIBILITY ;--------------------------------------------------------------- subarr = atan(subarr, /phase) * 180./!pi ; SIMPLE INTERPOLATION OF THE PHASES: PROBLEMATIC AT PHASE WRAPS ;--------------------------------------------------------------- ; result.visphi[i] = interpolate(subarr, xx, yy, cubic=-0.5) ; SELECT ELEMENTS WHICH SEEM TO HAVE WRAPPED AND REPLACE THEM ;--------------------------------------------------------------- select = where(abs(subarr-(subarr[0])) gt 200.) if select[0] ne -1 then begin subarr[select] += (subarr[0]/abs(subarr[0])*360.) endif ; DO INTERPOLATION ON THE NEWLY CREATED SUBARRAY ;--------------------------------------------------------------- result[j].visphi[i] = interpolate(subarr, xx, yy, cub=-0.5) endfor ; CALCULATE THE VISIBILITY ;----------------------------------------------------------------------- result.visamp[i] = result.corflx[i] / totflx ; CATCH RANDOM PHASES FOR EXTREMELY LOW VISIBILITIES ;----------------------------------------------------------------------- selidx = where(result.visamp[i] lt 1.0E-6) if selidx[0] ge 0 then result[selidx].visphi[i] = 0.0 ; COPY THE WAVELENGTH TO THE RESULT ;----------------------------------------------------------------------- result.wavele[i] = uvplan[i].wavele ; PLOT THE UV POINTS IN THE UV PLANE ;----------------------------------------------------------------------- if keyword_set(uvplot) then begin maxpbl = 2.0*max(result.pbline) modl_m2plot, uvplan[i], /fftd, flux=0, maxpbl=maxpbl, $ overpl=overpl if 0 then begin loadct, 1, /silent for j=0,n_elements(result)-1 do begin plots, result[j].ucoord, result[j].vcoord, psym=8, color=254*uvcoor[j].visamp[i] plots, -result[j].ucoord, -result[j].vcoord, psym=8, color=254*uvcoor[j].visamp[i] endfor endif else begin loadct, 39, /silent oplot, result.ucoord, result.vcoord, psym=8, color=50 oplot, -result.ucoord, -result.vcoord, psym=8, color=50 endelse overpl = 1 wait, 0.5 endif endfor ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, reform(result) END ;+ ; NAME: ; MIDI_DATCOM ; ; PURPOSE: ; Interpolate MIDI data at certain wavelengths. ; ; CATEGORY: ; MIDI model comparison. ; ; CALLING SEQUENCE: ; result = midi_datcom(caldat, wavele [, smtval]) ; ; REQUIRED INPUTS: ; caldat a MIDI result structure containing the original data ; wavele array specifying the wavelengths to be interpolated ; ; OPTIONAL INPUTS: ; smtval number of elements by which the data is to be smoothed ; ; KEYWORDS: ; wavwid half width of wavelength averaging window ; ; OUTPUTS: ; result structure containing the data at the selected wavelengths ; ; DESCRIPTION AND EXAMPLE: ; This function interpolates MIDI data at the wavelengths specified in ; WAVELE and returns the result in a similar structure. Before the ; interpolation is carried out, the data can be smoothed by setting - ; SMTVAL to a nonzero value. The program is invoked by: ; result = modl_datcom(caldat) ; ; CALLED BY: ; midi_2bbfit ; midi_3bbfit ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2009-07-14 Written by Konrad R. W. Tristram ; 2010-09-23 K. Tristram: Make result and input structs compatible ; 2010-12-02 K. Tristram: Make more flexible when copying from old tags ; 2012-10-17 K. Tristram: Major rewrite, added keyword WAVWID. ; FUNCTION MIDI_DATCOM, CALDAT, WAVELE, SMTVAL, WAVWID=WAVWID ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- if n_elements(wavele) lt 1 then wavele = 12.0 wavelm = n_elements(wavele) if n_elements(smtval) lt 1 then smtval = 0.0 if n_elements(wavwid) lt 1 then wavwid = 0.0 if n_elements(wavwid) ne wavelm then wwidth = replicate(wavwid, wavelm) ; GET THE INDICES OF THE DATA TAGS ;------------------------------------------------------------------------------- datnam = ['MSKFLX', 'MSKFLXERR', 'TOTFLX', 'TOTFLXERR', 'CORFLX', 'CORFLXERR', $ 'VISAMP', 'VISAMPERR', 'VISPHI', 'VISPHIERR'] tagnam = tag_names(caldat) datidx = intarr(n_elements(datnam)) for i=0,n_elements(datnam)-1 do datidx[i] = where(datnam[i] eq tagnam) datidx = datidx[where(datidx ge 0)] ;CREATE THE STRUCTURE HOLDING THE RESULT ;------------------------------------------------------------------------------- result = caldat[0] result = mod_struct(result, 'WAVELE', wavele) for i=0,n_elements(datidx)-1 do begin result = mod_struct(result, tagnam[datidx[i]], fltarr(wavelm)) endfor result = replicate(result, n_elements(caldat)) ; COPY BASELINE INFORMATION AND STRING VARIABLES ;------------------------------------------------------------------------------- for i=0,n_elements(tagnam)-1 do $ if (total(i eq datidx) ne 1) and tagnam[i] ne 'WAVELE' then begin result.(i) = caldat.(i) endif ; IF WAVELENGTH WINDOW IS ZERO, THEN HAVE TO INTERPOLATE ;------------------------------------------------------------------------------- if wavwid eq 0.0 then begin ; LOOP OVER DATA SETS AND OVER THE DATA TAGS ;------------------------------------------------------------------------------- for i=0,n_elements(caldat)-1 do for j=0,n_elements(datidx)-1 do begin result[i].(datidx[j]) = interpol(smooth(caldat[i].(datidx[j]), smtval, /NaN), $ caldat[i].wavele, wavele, /quadr) endfor ; WAVELENGTH WINDOW IS NOT ZERO ;------------------------------------------------------------------------------- endif else begin ; LOOP OVER DATA SETS AND WAVELENGTHS AND SELECT THE WAVELENGTH ELEMENTS ;------------------------------------------------------------------------------- for i=0,n_elements(caldat)-1 do for k=0,n_elements(wavele)-1 do begin wavsel = where((caldat[i].wavele gt (wavele[k]-wwidth[k])) and $ (caldat[i].wavele lt (wavele[k]+wwidth[k]))) ; LOOP OVER THE DATA TAGS ;----------------------------------------------------------------------- for j=0,n_elements(datidx)-1 do begin ; CALCULATE THE SMOOTHED SPECTRUM ;--------------------------------------------------------------- tmpvar = smooth(caldat[i].(datidx[j]), smtval, /NaN) ; IF THERE ARE SELECTED WAVELENGTH BINS THEN AVERAGE THEM ;--------------------------------------------------------------- if wavsel[0] ge 0 then begin result[i].(datidx[j])[k] = mean(tmpvar[wavsel]) ; OTHERWISE HAVE TO INTERPOLATE ;--------------------------------------------------------------- endif else begin result[i].(datidx[j])[k] = interpol(tmpvar, $ caldat[i].wavele, wavele[k], /quadr) endelse endfor endfor endelse ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, result END ;+ ; NAME: ; MIDI_DPHASE ; ; PURPOSE: ; Determine the differential phase from input phase data ; ; CATEGORY: ; MIDI model comparison. ; ; CALLING SEQUENCE: ; result = modl_dphase(mdlvis [, method] [, wavran=wavran] ; [, smtval=smtval] [, /output]) ; ; REQUIRED INPUTS: ; mdlvis a MIDI result structure containing the original data ; ; OPTIONAL INPUTS: ; method method to be used for determining the differential phase: ; 0: Group delay method using EWS ; 1: AMBER method ; 2: linear fit to the phase in k-space ; ; KEYWORDS: ; wavran two element array specifying the wavelength range ; smtval smooth the phases before doing any calculations ; output make plots comparing original to differential phases ; silent suppress output of most warning messages ; ; OUTPUTS: ; result structure containing the data with differential phases only ; ; DESCRIPTION AND EXAMPLE: ; This function calculates the differential phases from input phase ; data, that is the constant and linear components in the phases are ; removed. The input structure must have the form {wavele:wavele, ; visphi:visphi, visphierr:visphierr}, e.g. it must be a MIDI result ; structure, and it must contain at least 3 wavelength settings. There ; are three methods for determining the differential phase. Method 0 ; uses the algorithm implemented in EWS, where the additional delay ; is determined from the group delay. Method 1 uses the method by ; Florentin Millour for AMBER, where the slope in the phase is deter- ; mined by differentiation of the linear expansion of the phase. Doing ; the calculation in the complex plane removes problems due to the ; wrapping of the phase. This method has problems for noise on the ; order of !pi. In method 2, the phases are unwrapped be looking for ; jumps on the order of 2*!pi and then rejoining the segments. Then a ; simple linear fit is carried out. The wavelength range to be used for ; the determination of the constant and linear component can be speci- ; fied by the keyword WAVRAN; finally, the correction is applied to all ; wavelengths. With the keyword SMTVAL set, the phases are smoothed ; before some of the calculations are done. ; The programme is called by e.g.: result = modl_dphase(mdlvis) ; ; CALLED BY: ; midi_2bbdat ; midi_2bbcmp ; midi_3bbdat ; midi_3bbcmp ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2009-07-17 Written by Konrad R. W. Tristram ; 2009-07-20 K. Tristram: Added the group delay method from EWS. ; 2009-08-17 K. Tristram: Catch random phases for very low visibility. ; 2009-08-18 K. Tristram: Added keyword SILENT and set allphases to 0 ; for less than three wavelength elements. ; FUNCTION MIDI_DPHASE, MDLVIS, METHOD, WAVRAN=WAVRAN, SMTVAL=SMTVAL, $ OUTPUT=OUTPUT, SILENT=SILENT ; INITIALISATION OF SEVERAL VARIABLES ;------------------------------------------------------------------------------- if n_elements(wavran) ne 2 then wavran = [8.0, 13.2] else wavran = float(wavran) if n_elements(method) le 0 then method = 0 if n_elements(smtval) le 0 then smtval = 0 dummyv = 'string' imagin = complex(0,1) ; CHECK IF THERE ARE AT LEAST 3 WAVELENGTHS ELEMENTS IN THE INPUT STRUCTURE ;------------------------------------------------------------------------------- if n_elements(mdlvis[0].wavele) lt 3 then begin print, 'Not enough wavelength bins, need at least 3. Set phases to 0.' result = mdlvis result.visphi = 0.0 return, result endif ; CREATE RESULT STRUCTURE AS A COPY OF THE INPUT STRUCTURE ;------------------------------------------------------------------------------- result = mdlvis ; LOOP OVER THE UV POINTS IN THE RESULT ;------------------------------------------------------------------------------- for i=0,n_elements(result)-1 do begin ; GET THE WAVELENGTHS AND SELECT THE ONES WITHIN THE WAVELENGTH RANGE ;----------------------------------------------------------------------- wavele = result[i].wavele selidx = where((wavele ge wavran[0]) and (wavele le wavran[1])) ; CHECK IF THERE ARE ENOUGH ELEMENTS, OTHERWISE SKIP THIS UV POINT ;----------------------------------------------------------------------- if n_elements(selidx) lt 3 then begin if ~ keyword_set(silent) then begin print, 'Not enough wavelength bins, need at ' + $ 'least 3. Skipping UV point # ' + $ string(i, format='(I4)') + '.' endif endif else begin ; CHECK FOR VERY LOW VISIBILITIES WHICH HAVE USELESS PHASE VALUES ;----------------------------------------------------------------------- tmpval = where(result[i].visamp[selidx] gt 1.0E-6) if n_elements(tmpval) lt 3 then begin if ~ keyword_set(silent) then begin print, 'Less than 3 points with visibility ' + $ '> 1.0E-6. Skipping UV point # ' + $ string(i, format='(I4)') + '.' endif endif else begin ; DETERMINE THE FINAL SELCTION OF WAVELENGTH BINS ;----------------------------------------------------------------------- selidx = selidx[tmpval] if selidx[0] eq 0 then selidx = selidx[1:*] ; CALCULATE THE WAVE NUMBER AND GET THE INPUT PHASES ;----------------------------------------------------------------------- wavenu = 1.0 / wavele phaval = result[i].visphi phaerr = result[i].visphierr > 0.001 ; SELECT THE METHOD FOR DOING THE REMOVAL ;----------------------------------------------------------------------- case method of ; METHOD 1: USE THE GROUP DELAY METHOD FROM EWS ;--------------------------------------------------------------- 0: begin ; COMPUTE THE COMPLEX VISIBILITY AND THE RANGE OF K ;------------------------------------------------------- visold = result[i].corflx*exp(imagin*!pi/180.*phaval) k_valu = 2 * !pi * wavenu[selidx] maxi_k = max(k_valu) ; CREATE A FOURIER TRANSFORM OBJECT ;------------------------------------------------------- fftobj = obj_new('fft1d', [0, 20*maxi_k/4096.,4096]) ; GET THE GRID OF DELAYS AND THE DELAY FUNCTION ;------------------------------------------------------- delgri = fftobj->kgrid() delfun = fftobj->cfft(k_valu, visold[selidx]) ; DESTROY THE FOURIER TRANSFORM OBJECT ;------------------------------------------------------- obj_destroy, fftobj ; GET THE PEAK OF THE DELAY AND REMOVE THE DELAY ;------------------------------------------------------- delpea = pk(delgri, abs(delfun)) visnew = visold * exp(imagin * 2*!pi*wavenu * delpea) ; ALSO REMOVE THE PHASE OFFSET ;------------------------------------------------------- vistot = conj(total(visnew[selidx])) visnew = visnew * (vistot / abs(vistot)) phanew = atan(visnew, /phase) / !pi * 180 ; CALCULATE THE PHASE COMPONENT THAT WAS REMOVED ;------------------------------------------------------- phafit = exp(imagin * 2*!pi*wavenu * delpea) phafit = atan(phafit, /phase) / !pi * 180 end ; METHOD 1: USE THE METHOD USED FOR AMBER (FROM F. MILLOUR) ;--------------------------------------------------------------- 1: begin ; OPTIONALLY SMOOTH THE PHASES AND ERRORS ;------------------------------------------------------- phaval = smooth(phaval, smtval) phaerr = smooth(phaerr, smtval) ; GET MEAN DELAY OF THE PHASE IN THE COMPLEX PLANE ;------------------------------------------------------- tmpvar = exp(imagin*!pi/180.*(phaval-shift(phaval, 1))) tmpvar = atan(tmpvar, /phase)/(wavenu-shift(wavenu, 1)) delavg = mean(tmpvar[selidx]) ; REMOVE THIS DELAY FROM THE ORIGINAL PHASE ;------------------------------------------------------- phanew = !pi/180. * phaval - delavg * wavenu phanew = exp(imagin * (phanew - mean(phanew[selidx]))) phanew = atan(phanew, /phase) / !pi * 180 phanew = phanew - mean(phanew[selidx]) ; CALCULATE THE PHASE COMPONENT THAT WAS REMOVED ;------------------------------------------------------- phafit = exp(imagin * delavg * wavenu) phafit = atan(phafit, /phase) / !pi * 180 end ; METHOD 2: DO A LINEAR FIT TO THE PHASES AFTER UNWRAPPING THEM ;--------------------------------------------------------------- else: begin ; REMOVE PHASE WRAPPING BY ADDING +/-180° AT JUMPS ;------------------------------------------------------- differ = phaval - shift(phaval,1) jmpidx = (where((abs(differ) gt 260)[1:*]) + 1)[0] while jmpidx ge 1 do begin phaval[jmpidx:*] -= differ[jmpidx] / $ abs(differ[jmpidx]) * 360. differ = phaval - shift(phaval,1) jmpidx = (where((abs(differ) gt 260)[1:*])+1)[0] endwhile ; OPTIONALLY SMOOTH THE PHASES AND ERRORS ;------------------------------------------------------- phaval = smooth(phaval, smtval) phaerr = smooth(phaerr, smtval) ; MAKE A LINEAR FIT TO THE DATA AND REMOVE IT ;------------------------------------------------------- fitpar = linfit(wavenu[selidx], phaval[selidx], $ measure=phaerr[selidx]) phanew = phaval - fitpar[0] - fitpar[1] * wavenu ; CALCULATE THE PHASE COMPONENT REMOVED ;------------------------------------------------------- phafit = fitpar[0] + fitpar[1] * wavenu end endcase ; COPY THE NEW PHASES TO THE RESULT STRUCTURE ;----------------------------------------------------------------------- result[i].visphi = phanew endelse ; SET THE PHASES TO ZERO WHERE THE VISIBILITY IS EXTREMELY LOW ;----------------------------------------------------------------------- selidx = where(result[i].visamp le 1.0E-6) if selidx[0] ge 0 then result[i].visphi[selidx] = 0.0 ; OPTIONALLY MAKE A PLOT OF THE ORIGINAL DATA, THE FIT AND THE RESULT ;----------------------------------------------------------------------- if keyword_set(output) then begin plot, 1/reverse(wavran), [min((phaval-phaerr)[selidx]) < 0, $ max((phaval+phaerr)[selidx]) > 0], xstyle=1, /nodata, $ xtitle='wavenumber k = 1/!4k!X [1/!4l!Xm]', $ ytitle='differential phase [' + string("260B) + ']' oplot, [0,1], [0,0], linestyle=1 oplot, wavenu, phaval, psym=-2, color=90 errplot, wavenu, phaval-phaerr, phaval+phaerr, color=90, width=0 oplot, wavenu[selidx], phafit[selidx], color=150 oplot, wavenu, phanew, color=70 read, dummyv endif endelse endfor ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, result END ;+ ; NAME: ; MIDI_PLTCMP ; ; PURPOSE: ; Plot a comparison of measured data and modelled data. ; ; CATEGORY: ; MIDI model comparison. ; ; CALLING SEQUENCE: ; modl_pltcmp, datvis, mdlvis [, mdlvi1] [, mdlvi2] [, wavran=wavran] ; [, /corflx] [, /visamp] [, /visphi] [, /normal] ; [, pmulti=pmulti] [, psfile=psfile] [, /compact] ; ; REQUIRED INPUTS: ; datvis a MIDI result structure containing the measured data ; mdlvis a MIDI result structure containing the model data ; ; OPTIONAL INPUTS: ; mdlvi1 a MIDI result structure with data of a second model ; mdlvi2 a MIDI result structure with data of a third model ; ; KEYWORDS: ; wavran the wavelength range for which the plots will be produced ; corflx make comparison plots for the correlated fluxes ; visamp make comparison plots for the visibilities ; visphi make comparison plots for the differential phases ; normal normalise the scaling of the y axis in all plots ; pmulti number of colunms and rows per page in the plot ; ewidth width of the error bar ends ; psfile prefix of the PS file for POSTSCRIPT output ; compact make a compact plot in the case of PS output ; ; OUTPUTS: ; - plots comparing the measured data to the model data ; ; DESCRIPTION AND EXAMPLE: ; This programme creates plots of the correlated flux, the visibility ; and / or the phase as a function of wavelength, comparing measured ; data and model data. A comparison plot is produced for every UV ; point. Without that the keywords CORFLX, VISAMP and VISPHI are set, ; only plots for the correlated flux and the phase are created. Using ; the keywords MDLVI1 and MDLVI2 two additional models can be provided, ; for example the contributions of the two individual components of a ; two component model. The programme can be either used for on screen ; display of the plots or for creating POSTSCRIPT plots. The programme ; is simply called by: midi_pltcmp, datvis, mdlvis ; ; CALLED BY: ; midi_2bbcmp ; midi_3bbcmp ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2007-02-07 Adapted from earlier test_5.pro by Konrad R. Tristram ; 2007-10-22 K. Tristram: Changed y-axis scaling of plots. ; 2008-09-17 K. Tristram: Adaptation to the new MIDI result format. ; 2009-07-20 K. Tristram: Completely rewritten and programme renamed. ; 2009-08-18 K. Tristram: Added keyword PMULTI. ; 2012-11-19 K. Tristram: Added keyword EWIDTH and limit yrange of plots. ; PRO MIDI_PLTCMP, DATVIS, MDLVIS, MDLVI1, MDLVI2, WAVRAN=WAVRAN, CORFLX=CORFLX, $ VISAMP=VISAMP, VISPHI=VISPHI, NORMAL=NORMAL, PMULTI=PMULTI, $ EWIDTH=EWIDTH, PSFILE=PSFILE, COMPACT=COMPACT ; CHECK THE INPUT VARIABLES AND SET THEIR DEFAULT VALUES ;------------------------------------------------------------------------------- if n_elements(ewidth) le 0 then ewidth = 0.0 if n_elements(corflx) le 0 then corflx = 1 if n_elements(visamp) le 0 then visamp = 0 if n_elements(visphi) le 0 then visphi = 1 if n_elements(wavran) ne 2 then xrange = [7.85, 13.25] else xrange = wavran if keyword_set(psfile) and keyword_set(compact) then normal = 1 ; INITIALISE SEVERAL OTHER VARIABLES ;------------------------------------------------------------------------------- x_titl = 'wavelength !4k!X [!4l!Xm]' y_titl = ['correlated flux F!Dcor!N [Jy]', 'visibility V', 'phase !4u!X [' + $ string("260B) + ']'] msgcoo = 'Warning: UV points do not seem to match, euclidean distance in ' + $ 'UV plane more than 1 m.' dummyv = 'string' notick = replicate(' ', 20) numele = n_elements(mdlvis) do_plt = [keyword_set(corflx), keyword_set(visamp), keyword_set(visphi)] factor = 1.0 ; CHECK IF THE NUMBER OF UV POINTS MATCH, OTHERWISE EXIT ;------------------------------------------------------------------------------- if n_elements(datvis) ne numele then begin message, 'Number of UV points do not match. Returning.', /continue return endif ; CHECK IF ANY PLOTS ARE SELECTED AT ALL, OTHERWISE EXIT ;------------------------------------------------------------------------------- if total(do_plt) le 0 then begin message, 'No type of plot chosen. Returning.', /continue return endif ; CHECK IF THERE IS MORE THAN ONE MODEL TO PLOT ;------------------------------------------------------------------------------- plt_m2 = (plt_m1 = 0) if (n_elements(mdlvi1) eq numele) and (size(mdlvi1, /type) eq 8) then plt_m1 = 1 if (n_elements(mdlvi2) eq numele) and (size(mdlvi2, /type) eq 8) then plt_m2 = 1 ; SELECT THE INDICES OF ELEMENTS WITHIN THE WAVELENGTH RANGE ;------------------------------------------------------------------------------- ds = where((datvis.wavele gt xrange[0]) and (datvis.wavele lt xrange[1])) ms = where((mdlvis.wavele gt xrange[0]) and (mdlvis.wavele lt xrange[1])) ; CHECK IF THERE ARE ENOUGH ELEMENTS WITHIN THE WAVELENGTH RANGE ;------------------------------------------------------------------------------- if (ds[0] lt 0) or (ms[0] lt 0) then begin message, 'Not enough elements in the wavelength range. Returning.', $ /continue return endif ; CREATE THE ARRAYS FOR THE DATA TO BE PLOTTED ;------------------------------------------------------------------------------- tmpvar = !values.f_nan data_y = make_array(n_elements(datvis[0].wavele), numele, 3, value=tmpvar) data_e = data_y mdl0_y = make_array(n_elements(mdlvis[0].wavele), numele, 3, value=tmpvar) if plt_m1 then mdl1_y = make_array(n_elements(mdlvi1[0].wavele), numele, $ 3, value=tmpvar) $ else mdl1_y = make_array(3, numele, 3, value=tmpvar) if plt_m2 then mdl2_y = make_array(n_elements(mdlvi2[0].wavele), numele, $ 3, value=tmpvar) $ else mdl2_y = make_array(3, numele, 3, value=tmpvar) yrange = make_array(2, 3, value=0.0) ; DEPENDING ON WHICH PLOTS ARE SELECTED, COPY THE DATA INTO THESE ARRAYS ;------------------------------------------------------------------------------- if corflx then begin data_y[*,*,0] = datvis.corflx data_e[*,*,0] = datvis.corflxerr mdl0_y[*,*,0] = mdlvis.corflx if plt_m1 then mdl1_y[*,*,0] = mdlvi1.corflx if plt_m2 then mdl2_y[*,*,0] = mdlvi2.corflx endif if visamp then begin data_y[*,*,1] = datvis.visamp data_e[*,*,1] = datvis.visamperr mdl0_y[*,*,1] = mdlvis.visamp if plt_m1 then mdl1_y[*,*,1] = mdlvi1.visamp if plt_m2 then mdl2_y[*,*,1] = mdlvi2.visamp endif if visphi then begin data_y[*,*,2] = datvis.visphi data_e[*,*,2] = datvis.visphierr mdl0_y[*,*,2] = mdlvis.visphi if plt_m1 then mdl1_y[*,*,2] = mdlvi1.visphi if plt_m2 then mdl2_y[*,*,2] = mdlvi2.visphi endif if plt_m1 then mdl1_x = mdlvi1.wavele $ else mdl1_x = make_array(3, numele, value=tmpvar) if plt_m2 then mdl2_x = mdlvi1.wavele $ else mdl2_x = make_array(3, numele, value=tmpvar) ; IF PLOTS ARE TO BE NORMALISED, THEN LOOP OVER THE TYPE OF PLOTS ;------------------------------------------------------------------------------- if keyword_set(normal) then for j=0,2 do if do_plt[j] then begin ; CALCULATE GENERAL YRANGES FOR THE PLOTS ;----------------------------------------------------------------------- yrange[*,j] = [min(smooth([data_y[ds,*,j]-data_e[ds,*,j], $ mdl0_y[ms,*,j]], 3)), $ max(smooth([data_y[ds,*,j]+data_e[ds,*,j], $ mdl0_y[ms,*,j]], 3))] ; ADD A 5 PERCENT MARGIN TO THE DATA RANGE ;----------------------------------------------------------------------- yrange[0,j] -= 0.50 * abs(yrange[0,j]) yrange[1,j] += 0.50 * abs(yrange[1,j]) ; FOR THE VISIBILITY AND PHASES LIMIT THE DATA RANGE ;----------------------------------------------------------------------- if j eq 0 then yrange[*,j] = yrange[*,j] > 0 if j eq 1 then yrange[*,j] = yrange[*,j] > 0 < 2 if j eq 2 then yrange[*,j] = yrange[*,j] > (-190) < 190 endif ; SAVE THE ORIGINAL PLOT SETTINGS AND LOAD COLOR TABLE 39 ;------------------------------------------------------------------------------- d_name = !d.name p = !p omargi = [!x.omargin, !y.omargin] loadct, 39, /silent ; ENTER ON SCREEN PLOT SECTION IF NO PS FILENAME IS GIVEN ;------------------------------------------------------------------------------- if ~ keyword_set(psfile) then begin ; SET THE PLOT DEVICE TO THE SCREEN ;----------------------------------------------------------------------- if !version.os_family eq 'Windows' then set_plot, 'WIN' $ else set_plot, 'X' device, decompose = 0 ; SET THE MARGINS FOR THE PLOTS AND THE CHARACTER SIZE ;----------------------------------------------------------------------- x_marg = [8.0,1.0] y_marg = [3.2,2.0] !x.omargin = [0.0, 0.0] !y.omargin = [0.0, 0.0] !p.charsize = 1.25 ; SET THE MULTIPLOT ENVIRONMENT DEPENDING ON THE NUMBER OF PLOTS ;----------------------------------------------------------------------- case total(do_plt) of 1: !p.multi = [0,1,1] 2: !p.multi = [0,2,1] else: begin !p.multi = [0,3,1] !p.charsize = 2. * !p.charsize end endcase ; OPEN A NEW PLOT WINDOW ;----------------------------------------------------------------------- winsiz = factor * !p.multi[1:2] * [350, 350] window, xsize=winsiz[0], ysize=winsiz[1], /free ; LOOP OVER THE UV POINTS ;----------------------------------------------------------------------- for i=0,numele-1 do begin ; CALCULATE THE EUCLIDEAN DISTANCE OF MODEL & DATA UV COORDS ;--------------------------------------------------------------- uvdist = sqrt((datvis[i].ucoord - mdlvis[i].ucoord)^2 + $ (datvis[i].ucoord - mdlvis[i].ucoord)^2) ; PRINT OUT A WARNING IF THE COORDINATES ARE TOO DIFFERENT ;--------------------------------------------------------------- if uvdist gt 1 then print, msgcoo ; SELECT THE INDICES OF ELEMENTS WITHIN THE WAVELENGTH RANGE ;--------------------------------------------------------------- ds = where((datvis[i].wavele gt xrange[0]) and $ (datvis[i].wavele lt xrange[1])) ms = where((mdlvis[i].wavele gt xrange[0]) and $ (mdlvis[i].wavele lt xrange[1])) ; CHECK IF THERE ARE ENOUGH ELEMENTS WITHIN THE WAVELENGTH RANGE ;--------------------------------------------------------------- if (ds[0] lt 0) or (ms[0] lt 0) then begin message, 'Not enough elements in the wavelength ra' + $ 'nge. Skipping UV point # '+strcompress(i)+'.', /cont endif else begin ; LOOP OVER THE THREE POSSIBLE PLOTS AND CHECK WHICH ONES TO DO ;--------------------------------------------------------------- for j=0,2 do if do_plt[j] then begin ; CALCULATE THE DATA RANGE ;------------------------------------------------------- datran = [min(smooth([data_y[ds,i,j]-data_e[ds,i,j], $ mdl0_y[ms,i,j]], 3)), $ max(smooth([data_y[ds,i,j]+data_e[ds,i,j], $ mdl0_y[ms,i,j]], 3))] ; ADD A 5 PERCENT MARGIN TO THE DATA RANGE ;------------------------------------------------------- datran[0] -= 0.1 * abs(datran[0]) datran[1] += 0.1 * abs(datran[1]) ; FOR THE VISIBILITY AND PHASES LIMIT THE DATA RANGE ;------------------------------------------------------- if j eq 1 then datran = datran > 0 < 2 if j eq 2 then datran = datran > (-190) < 190 ; CREATE THE PLOT AXIS AND AXIS LABELS ;------------------------------------------------------- plot, xrange, datran, xrange=xrange, xstyle=1, $ yrange=yrange[*,j], ystyle=ystyle, ytitle=y_titl[j], $ xtitle=x_titl, xmargin=x_marg, ymargin=y_marg, /nodata ; FOR VISIBILITY AND PHASE PLOT HELPER LINES ;------------------------------------------------------- if j eq 1 then oplot, [0,100], [1,1], linestyle=1 if j eq 2 then oplot, [0,100], [0,0], linestyle=1 ; PLOT THE DATA WITH ERROR BARS AND THE MODELS ;------------------------------------------------------- oplot, datvis[i].wavele, data_y[*,i,j], color=150 errplot, datvis[i].wavele, width=ewidth, $ data_y[*,i,j]-data_e[*,i,j], $ data_y[*,i,j]+data_e[*,i,j], color=150 oplot, mdlvis[i].wavele, mdl0_y[*,i,j], color=60 oplot, mdl1_x[*,i], mdl1_y[*,i,j], color=60, linestyle=1 oplot, mdl2_x[*,i], mdl2_y[*,i,j], color=60, linestyle=2 endif ; CREATE THE TITLE STRING FOR THE PLOTS ;--------------------------------------------------------------- titstr = 'UV point'+string(i, format='(I4)')+': u =' + $ string(datvis[i].ucoord, format='(F6.1)')+' m, v ='+$ string(datvis[i].vcoord, format='(F6.1)')+' m' ; PUT MORE INFORMATION IN THE TITLE IF THERE ARE MORE PANELS ;--------------------------------------------------------------- if total(do_plt) gt 1 then titstr += ', BL =' + $ string(datvis[i].pbline, format='(F6.1)')+' m, PA ='+ $ string(datvis[i].pangle, format='(F6.1)') + $ string("260B) + ' (' + datvis[i].baseli + ')' ; PLOT THE TITLE ;--------------------------------------------------------------- xyouts, winsiz[0]/2, winsiz[1]-15, titstr, align=0.5, $ charsize=1.5, /device ; READ A CHARACTER BEFORE DOING THE NEXT PLOT OR EXITING ;--------------------------------------------------------------- read, 'Press a Enter to continue or "q" to quit.', dummyv if total(strupcase(dummyv) eq ['F', 'Q', 'QUIT']) then begin set_plot, d_name !p = p return endif endelse ; if (ds[0] lt 0) or (ms[0] lt 0) then else begin endfor ; ENTER OUTPUT TO PS FILE IF A PS FILENAME IS GIVEN ;------------------------------------------------------------------------------- endif else begin ; SET THE PLOT DEVICE TO THE PS DEVICE ;----------------------------------------------------------------------- set_plot, 'ps' ; DEFINE THE SUFFIXES OF THE PS FILES ;----------------------------------------------------------------------- ending = ['cor.ps', 'vis.ps', 'phi.ps'] ; SET THE MULTIPLOT ENVIRONMENT ;----------------------------------------------------------------------- if n_elements(pmulti) ne 2 then !p.multi = [0, 4, 6] $ else !p.multi = [0, pmulti] numplt = !p.multi[1] * !p.multi[2] ; SET THE MARGINS FOR THE PLOTS AND THE CHARACTER SIZE ;----------------------------------------------------------------------- charsz1 = 1.5 if keyword_set(compact) then begin xmargin = [0.0, 0.0] ymargin = [0.0, 0.0] charsz2 = 2.5/!p.multi[1] !x.omargin = [8.0, 8.0] !y.omargin = [5.0, 1.0] endif else begin xmargin = [8.0, 1.0] ymargin = [3.2, 1.0] charsz2 = 2.0/!p.multi[1] !x.omargin = [0.0, 0.0] !y.omargin = [2.0, 0.0] endelse ; LOOP OVER THE THREE POSSIBLE PLOTS AND CHECK WHICH ONES TO DO ;----------------------------------------------------------------------- for j=0,2 do if do_plt[j] then begin ; START THE MULTIPLE PLOTS AT THE TOP OF THE PAGE ;----------------------------------------------------------------------- !p.multi[0] = 0 ; OPEN A NEW PS FILE ;----------------------------------------------------------------------- device, xsize=19.0, ysize=27.0, xoffset=1.0, yoffset=1.3, $ file=psfile+ending[j], /portrait, /color ; LOOP OVER THE UV POINTS ;----------------------------------------------------------------------- for i=0,numele-1 do begin ; SELECT THE INDICES OF ELEMENTS IN THE WAVELENGTH RANGE ;--------------------------------------------------------------- ds = where((datvis[i].wavele gt xrange[0]) and $ (datvis[i].wavele lt xrange[1])) ms = where((mdlvis[i].wavele gt xrange[0]) and $ (mdlvis[i].wavele lt xrange[1])) ; CHECK IF THERE ARE ENOUGH ELEMENTS IN THE RANGE ;--------------------------------------------------------------- if (ds[0] lt 0) or (ms[0] lt 0) then begin message, 'Not enough elements in the wave' + $ 'length range. Skipping UV point # '+ $ strcompress(i)+'.', /continue endif else begin ; CALCULATE THE DATA RANGE ;--------------------------------------------------------------- datran = [min(smooth([data_y[ds,i,j]-data_e[ds,i,j], $ mdl0_y[ms,i,j]], 3)), $ max(smooth([data_y[ds,i,j]+data_e[ds,i,j], $ mdl0_y[ms,i,j]], 3)) * 1.1] ; FOR THE VISIBILITY AND PHASES LIMIT THE DATA RANGE ;--------------------------------------------------------------- if j eq 0 then datran = datran > 0 if j eq 1 then datran = datran > 0 < 2 if j eq 2 then datran = datran > (-190) < 190 ; SET THE DEFAULT VALUES FOR THE AXIS TITLES AND TICK NAMES ;--------------------------------------------------------------- pltaxi = 0 xtitle = '' ytitle = '' ytickn = notick xtickn = notick ; DEPENDING ON WHERE PLOT IS ALTER THE TITLES AND TICK NAMES ;--------------------------------------------------------------- if keyword_set(compact) then case 1 of ; TOP LEFT ((!p.multi[0] eq 0) and $ ((numele - 1 - i) gt !p.multi[1])): begin ytitle = y_titl[j] ytickn = '' end ; TOP RIGHT (!p.multi[0] eq (!p.multi[1]*(!p.multi[2]-1)+1) and $ ((numele - 1 - i) ge !p.multi[1])): pltaxi = 1 ; BOTTOM LEFT ((!p.multi[0] eq !p.multi[1]) or $ (((!p.multi[0] mod !p.multi[1]) eq 0) and $ ((numele - 1 - i) lt !p.multi[1]))): begin xtickn = '' ytickn = '' xtitle = x_titl ytitle = y_titl[j] end ; BOTTOM RIGHT (((numele-1-i) eq ((numele-1 mod !p.multi[1])+1)) or $ (!p.multi[0] eq 1) or (i eq numele - 1)): begin xtitle = x_titl xtickn = '' pltaxi = 1 end ; BOTTOM MIDDLE (((!p.multi[0] gt 1) and (!p.multi[0] lt $ !p.multi[1])) or ((numele-1-i) lt !p.multi[1])): begin xtitle = x_titl xtickn = '' end ; LEFT MIDDLE (((!p.multi[0] mod !p.multi[1]) eq 0) and $ (!p.multi[0] gt !p.multi[1])): begin ytitle = y_titl[j] ytickn = '' end ; RIGHT MIDDLE ((((!p.multi[0]-1) mod !p.multi[1]) eq 0) and $ (!p.multi[0] lt (!p.multi[1]*(!p.multi[2]-1)+1)) $ and (!p.multi[0] gt 1)): pltaxi = 1 ; REST else: endcase else begin xtickn = '' ytickn = '' xtitle = x_titl ytitle = y_titl[j] endelse ; CREATE THE PLOT AXIS AND AXIS LABELS ;--------------------------------------------------------------- plot, xrange, datran, xrange=xrange, xstyle=1, $ yrange=yrange[*,j], xtitle=xtitle, ytitle=ytitle, $ xtickname=xtickn, ytickname=ytickn, xmargin=xmargin, $ ymargin=ymargin, charsize=charsz1, /nodata ; FOR VISIBILITY AND PHASE PLOT HELPER LINES ;--------------------------------------------------------------- if pltaxi then axis, yaxis=1, ystyle=1, ytitle=y_titl[j], $ charsize=charsz1 ; FOR VISIBILITY AND PHASE PLOT HELPER LINES ;--------------------------------------------------------------- if j eq 1 then oplot, [0,100], [1,1], linestyle=1 if j eq 2 then oplot, [0,100], [0,0], linestyle=1 ; PLOT THE DATA WITH ERROR BARS AND THE MODELS ;--------------------------------------------------------------- oplot, datvis[i].wavele, data_y[*,i,j], color=150 errplot, datvis[i].wavele, data_y[*,i,j]-data_e[*,i,j], $ data_y[*,i,j]+data_e[*,i,j], color=150, width=ewidth oplot, mdlvis[i].wavele, mdl0_y[*,i,j], color=60 oplot, mdl1_x[*,i], mdl1_y[*,i,j], color=60, linestyle=1 oplot, mdl2_x[*,i], mdl2_y[*,i,j], color=60, linestyle=2 ; GET THE TRUE PLOT RANGES ;--------------------------------------------------------------- x_min = !x.crange[0] x_ran = !x.crange[1]-!x.crange[0] y_min = !y.crange[0] y_ran = !y.crange[1]-!y.crange[0] ; PLOT SOME ADDITIONAL INFORMATION INTO THE PANEL ;--------------------------------------------------------------- xyouts, x_min + 0.05*x_ran, y_min + 0.93*y_ran, $ datvis[i].baseli, charsize=charsz2 xyouts, x_min + 0.95*x_ran, y_min + 0.93*y_ran, 'u =' + $ string(datvis[i].ucoord, format='(F5.1)') + ' m', $ align=1.0, charsize=charsz2 xyouts, x_min + 0.95*x_ran, y_min + 0.86*y_ran, 'v =' + $ string(datvis[i].vcoord, format='(F5.1)') + ' m', $ align=1.0, charsize=charsz2 ; PRINT OUT THE PAGE NUMBER ;--------------------------------------------------------------- if ((i mod (!p.multi[1]*!p.multi[2])) eq 0) then begin xyouts, 9500, 100, align=0.5, 'Page ' + $ string((i+1)/numplt+1, format='(I1)')+' of '+ $ string((numele-1)/numplt+1, format='(I1)'), $ charsize=1.0, /device endif endelse ; if (ds[0] lt 0) or (ms[0] lt 0) ... else begin endfor ; for i=0,numele-1 do begin ; CLOSE THE CURRENT PS FILE ;----------------------------------------------------------------------- device, /close endif ; if do_plt[j] then begin endelse ; RESTORE THE ORIGINAL PLOT SETTINGS ;------------------------------------------------------------------------------- set_plot, d_name !p = p !x.omargin = omargi[0:1] !y.omargin = omargi[2:3] END