; MIDI_DATSET__DEFINE - Automatically define a structure of type MIDI_DATSET. ; datset = {midi_datset} ; ; MIDI_GETKEY - Get FITS header keywords for a list of MIDI files. ; result = midi_getkey(filist, keywrd, direct=direct) ; ; MIDI_FILLIN - Automatically fill in information in a MIDI data set. ; midi_fillin, datset, /rmpath, /identi, /redpar, /caldat ; ; MIDI_GETDAT - Read raw and calibrated MIDI data from EWS into a struct. ; result = midi_getdat(datset, /rawdat, /fixbli, /silent) ; ; MIDI_CALINF - Get calibrator information for a list of MIDI results. ; calinf = midi_calinf(caldat, datset) ; ; MIDI_RVBCAL - Recalibrate fluxes with templates by Roy van Boekel. ; result = midi_rvbcal(caldat, calnam, calflx, /output) ; result = midi_rvbcal(caldat, calinf.calnam, calinf.calflx) ; ; MIDI_REDPLT - Plot the fluxes and visibilities of a MIDI result structure. ; midi_redplt, result, ctable, colour, /rawdat, psfile=psfile ; ; MIDI_MODPHT - Modify the sky bands and recompile EWS photometry programmes. ; midi_modpht, skypos, (/verbose) ; ; MIDI_REDSET - Fully reduce the data in a MIDI data set. ; midi_redset, prefix, datset, /noredu , /noplot, /shimsk, /slopef ; ; MIDI_PLTPRP - Plot a FITS header keyword for several MIDI data files. ; midi_pltprp, fileli, keywrd, ytitle, psfile=psfile ; ; MIDI_PLTPDS - Plot common keywords for fringe tracks in a MIDI data set. ; midi_pltpds, datset, titles, pspath=pspath ; ; MIDI_A2JAVE - Calculate the average and scatter of conversion factors. ; a2jave = midi_a2jave(adu2jy, mjdobs, nummax, sigclp=2.0) ; ; MIDI_ADU2JY - Calculate the ADU to JY conversion factors for calibrators. ; adu2jy = midi_adu2jy(datset, titles, pspath=pspath, /noplot) ; ; MIDI_CALERR - Add the errors from the variation of the transfer function. ; newdat = midi_calerr(caldat, datset, mincal=3, maxcal=5) ; ; MIDI_STRC2T - Print out an IDL structure as an ASCII table. ; midi_strc2t, datset, tagnam, maxwid=maxwid ; ; MIDI_DATA2T - Print out MIDI data as an ASCII file. ; midi_data2t, caldat, filnam, wavele ; ; MIDI_DATMRG - Merge several data sets using Gaussian error propagation ; result = midi_datmrg(caldat) ; ; MIDI_PLOT3D - Make a 3D plot of the UV plane using iPlot. ; midi_plot3d, caldat, wavele, ranges ;+ ; NAME: ; MIDI_DATSET__DEFINE ; ; PURPOSE: ; Automatically define a structure of type MIDI_DATSET. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; result = {midi_datset} ; ; REQUIRED INPUTS: ; None. ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORDS: ; None. ; ; OUTPUTS: ; result a structure of type MIDI_DATSET ; ; DESCRIPTION AND EXAMPLE: ; This procedure is for the automatic definition of a structure of type ; MIDI_DATSET. First the structure MIDI_REDPAR for the reduction ; parameters is defined, than the structure MIDI_DATSET itself is ; defined. The struct may then be created simply by calling: ; result = {midi_datset} ; ; CALLED BY: ; None. ; ; CALLING: ; None. ; ; MODIFICATION HISTORY: ; 2008-08-26 Written by Konrad R. Tristram ; 2010-02-05 Konrad R. Tristram: Added TAG for sky mask in REDPAR: SKYFIL. ; PRO MIDI_DATSET__DEFINE ; FIRST DEFINE THE REDUCTION PARAMETER STRUCT ;------------------------------------------------------------------------------- redpar = {midi_redpar, mskfil:'', skyfil:'', tsmoot:10.0, gsmoot:20.0, $ davera:1B, skydis:0.0, skypos:[0.0,0.0], nophot:1B, $ calibr:''} ; NOW DEFINE THE DATA SET STRUCT ITSELF ;------------------------------------------------------------------------------- datset = {midi_datset, datdir:'', identi:'', filnam:replicate('', 3), $ ewstag:'', dprcat:'', object:'', redpar:redpar, $ commen:'', calflx:0.0, caldia:0.0} END ;+ ; NAME: ; MIDI_GETKEY ; ; PURPOSE: ; Get FITS header keywords for a list of MIDI files. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; result = midi_getkey(filist, keywrd[, direct=direct]) ; ; REQUIRED INPUTS: ; filist a list of MIDI file names (individual or grouped) ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORDS: ; direct string specifiyng a directory where the files are located ; ; OUTPUTS: ; result table containing the file name used and the keyword values ; ; DESCRIPTION AND EXAMPLE: ; This function gets FITS header keywords for a list of MIDI files, which ; can be either a list of individual file names or a list of grouped file ; names. An optionally provided directory path is checked for the correct ; path separator and prefixed to the file name. The header keywords are ; then read using the IDL fits binary tables routines provided with ; MIA+EWS. The function requires the IDL fits binary tables routines to ; work. ; To get the object names for the fringe tracks in MIDI data sets, invoke ; result = midi_getkey(datset.filnam[0], 'object', direct=datset.datdir) ; ; CALLED BY: ; midi_fillin, midi_getdat, midi_pltres ; ; CALLING: ; None. ; ; MODIFICATION HISTORY: ; 2008-08-26 Written by Konrad R. Tristram ; FUNCTION MIDI_GETKEY, FILIST, KEYWRD, DIRECT=DIRECT ; CHECK IF DIRECTORY HAS BEEN SPECIFIED AND HOW MANY ELEMENTS ARE IN IT ;------------------------------------------------------------------------------- if (n_elements(direct) eq 1) then begin direct = replicate(string(direct), n_elements(filist)) endif if (n_elements(direct) ne n_elements(filist)) then begin direct = replicate('', n_elements(filist)) endif ; CREATE AN ARRAY HOLDING THE CORRECT FILE NAMES ;------------------------------------------------------------------------------- filnam = make_array(n_elements(filist), /string) ; LOOP OVER ALL ELEMENTS IN THE FILE LIST ;------------------------------------------------------------------------------- for i=0, n_elements(filist)-1 do begin ; REPLACE OR ADD THE CORRECT FILE SEPARATOR AT THE END OF THE DIRECTORY ;----------------------------------------------------------------------- lastch = strmid(direct[i], 0, 1, /reverse) if (lastch eq '\') || (lastch eq '/') then begin direct[i] = strmid(direct[i], 0, strlen(direct[i])-1) + $ path_sep() endif else if (direct[i] ne '') then direct[i] = direct[i] + path_sep() ; ONLY SELECT THE FIRST FILE IF THERE ARE SEVERAL IN COMPRESSED MODE ;----------------------------------------------------------------------- filnam[i] = direct[i] + (strsplit(filist[i], /extract))[0] ;print, filnam[i] endfor ; CREATE A FILE LIST STRUCTURE FROM THE FILES IN THE LIST ;------------------------------------------------------------------------------- filstr = obj_new('fitsfilelist', filnam) ; GET THE HEADER INFORMATION FROM THE FRINGE TRACK FILES ;------------------------------------------------------------------------------- filtbl = filstr->getheaderkeys(keywrd) ; DESTROY THE FILE OBJECTS ;------------------------------------------------------------------------------- obj_destroy, filstr ;RETURN THE RESULT ;------------------------------------------------------------------------------- return, filtbl END ;+ ; NAME: ; MIDI_FILLIN ; ; PURPOSE: ; Automatically fill in information in a MIDI data set. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; midi_fillin, datset [, /rmpath] [, /identi] [, /redpar] [, /caldat] ; ; REQUIRED INPUTS: ; datset a MIDI data set or list of MIDI data sets ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORDS: ; rmpath remove any previous path string in the file names ; identi also create an identity tag for the data sets ; redpar fill in default values for the reduction parameters ; caldat automatically add information from the calibrator database ; ; OUTPUTS: ; None. ; ; DESCRIPTION AND EXAMPLE: ; The programme completes the information in the provided data set(s). ; First, the directory path is added to the filenames if these are not ; present already. Then the filenames are used to extract FITS header ; keywords and to store them into the data set. Finally standard EWS ; tags are created. The programme is called by: ; midi_fillin, datset ; ; CALLED BY: ; None. ; ; CALLING: ; midi_getkey ; ; MODIFICATION HISTORY: ; 2008-08-26 Written by Konrad R. Tristram ; 2009-05-12 K. Tristram: Added keywords RMPATH, IDENTI & REDPAR. ; 2010-04-14 K. Tristram: Remove white spaces in object name. ; 2011-07-07 K. Tristram: Added search for calibrator information. ; 2011-08-18 K. Tristram: Added automatic calibrator matching. ; 2012-02-07 K. Tristram: Added handling of incomplete data sets. ; 2012-02-07 K. Tristram: Fixed fluxes taken from calibrator database. ; PRO MIDI_FILLIN, DATSET, RMPATH=RMPATH, IDENTI=IDENTI, REDPAR=REDPAR, $ CALDAT=CALDAT ; INITIALISE VARIABLES ;------------------------------------------------------------------------------- curdir = 'a/dummy\dir' if n_elements(rmpath) eq 1 then if rmpath eq 1 then rmpath = 'MIDI.' if n_elements(caldat) eq 0 then caldat = 1B ; SORT THE DATA SET ;------------------------------------------------------------------------------- datset = datset[sort(datset.datdir)] ; LOOP OVER THE ELEMENTS TO PREFIX THE DIRECTORY NAME TO THE FILE NAMES ;------------------------------------------------------------------------------- for i=0,n_elements(datset)-1 do begin ; REPLACE OR ADD THE CORRECT FILE SEPARATOR AT THE END OF THE DIRECTORY ;----------------------------------------------------------------------- lastch = strmid(datset[i].datdir, 0, 1, /reverse) if (lastch eq '\') || (lastch eq '/') then begin datset[i].datdir = strmid(datset[i].datdir, 0, $ strlen(datset[i].datdir)-1) + path_sep() endif else if (datset[i].datdir ne '') then begin datset[i].datdir = datset[i].datdir + path_sep() endif ; LOOP OVER THE THREE ELEMENTS OF THE FILENAME ARRAY ;----------------------------------------------------------------------- for j=0,2 do begin ; SPLIT STRING INTO ARRAY OF FILE NAMES ;--------------------------------------------------------------- tmpvar = strsplit(datset[i].filnam[j], /extr) ; LOOP OVER THE FILES IN THE ARRAY ;--------------------------------------------------------------- if (tmpvar[0] ne '') then for k=0,n_elements(tmpvar)-1 do begin ; REMOVE PATH ALREADY PRESENT IN FILE NAMES ;------------------------------------------------------- if keyword_set(rmpath) then begin nampos = strpos(tmpvar[k], rmpath) tmpvar[k] = strmid(tmpvar[k], nampos) endif ; CHECK WHETHER PATH IS ALRADY PREFIXED ;------------------------------------------------------- datdir = '*' + datset[i].datdir + '*' if not strmatch(tmpvar[k], datdir) then begin ; PREFIX PATH TO THE FILE NAME ;------------------------------------------------------- tmpvar[k] = datset[i].datdir + tmpvar[k] endif endfor datset[i].filnam[j] = strjoin(tmpvar, ' ') endfor endfor ; SET THE KEYWORDS THAT NEED TO BE READ ;------------------------------------------------------------------------------- keywrd = ['MJD-OBS', 'NRTS MODE', 'DPR CATG', 'DPR TYPE', 'OBS TARG NAME'] ; CONSTRUCT A FILE LIST FOR READING HEADER KEYWORDS ;------------------------------------------------------------------------------- filist = datset.filnam[0] for i=1,2 do begin tmpvar = where(filist eq '') if tmpvar[0] ge 0 then filist[tmpvar] = datset[tmpvar].filnam[i] endfor ; MAKE A TABLE CONTAINING THE HEADER INFORMATION ;------------------------------------------------------------------------------- filtbl = midi_getkey(filist, keywrd) ; COPY THE HEADER INFORMATION INTO THE STRUCT ;------------------------------------------------------------------------------- datset.dprcat = filtbl.dprcatg datset.object = strtrim(filtbl.obstargname, 2) ; LOAD THE CALIBRATOR DATA BASE OF ROY VAN BOEKEL ;------------------------------------------------------------------------------- if keyword_set(caldat) then boecal = vboekelbase() ; LOOP OVER THE ELEMENTS TO GENERATE THE TAGS AND GET THE CALIBRATOR DATA ;------------------------------------------------------------------------------- for i=0,n_elements(datset)-1 do begin ; FIND CALIBRATOR INFORMATION IF KEYWORD IS SET ;----------------------------------------------------------------------- if keyword_set(caldat) then begin ; FIND CALIBRATOR WITHIN 1 ARCMIN FROM POSITION IN FITS FILE ;--------------------------------------------------------------- srcdat = boecal->sourceFromFile(filist[i], rad=1./60., $ photData=phtdat, diamdata=diamet, ierr=errvar) ; IF A CALIBRATOR WAS FOUND FILL IN THE RESPECTIVE VALUES ;--------------------------------------------------------------- if not errvar then begin selbnd = [where(strmatch(phtdat.band, '*IRAS12*')), $ where(strmatch(phtdat.band, '*MSXC*'))] ; FIRST TRY THE IRAS 10µm FLUX: "IRAS12" ;------------------------------------------------------- if selbnd[0] ge 0 then begin datset[i].calflx = phtdat[selbnd[0]].flux ; THEN TRY THE MSX 12.13µm FLUX: "MSXC" ;------------------------------------------------------- endif else if selbnd[1] ge 0 then begin datset[i].calflx = phtdat[selbnd[1]].flux*1.44 ; OTHERWISE PUT A VERY NEGATIVE FLUX ;------------------------------------------------------- endif else datset[i].calflx = -999 ; FILL IN THE DIAMETER AND DPR CATEGORY ;------------------------------------------------------- datset[i].caldia = 1000.*diamet[0].diameter datset[i].dprcat = 'CALIB ' ; IF NO CALIBRATOR WAS FOUND THEN SET VALUES TO NAN ;--------------------------------------------------------------- endif else begin datset[i].calflx = !values.f_nan datset[i].caldia = !values.f_nan datset[i].dprcat = 'SCIENCE ' endelse endif ; IF THIS IS A NEW DIRECTORY THEN SET COUNTER J=0 ;----------------------------------------------------------------------- if curdir ne datset[i].datdir then begin curdir = datset[i].datdir j = 0 endif ; CREATE FIRST PART OF THE TAG DEPENDING ON THE DATA CATEGORY ;----------------------------------------------------------------------- if datset[i].dprcat eq 'SCIENCE ' then prefix = 'sci' else $ if datset[i].dprcat eq 'CALIB ' then prefix = 'cal' else begin message, 'Warning: DPR TYPE (SCIENCE/CALIB) not defined!', /cont prefix = 'unk' endelse ; ASSEMBLE THE FULL TAG NAME AND STORE IT INTO THE DATA SET ;----------------------------------------------------------------------- datset[i].ewstag = datset[i].datdir + 'ews/' + prefix + $ string(j, format='(I02)') ; ADD AN IDENTITY TAG FOR EVERY DATASET ;----------------------------------------------------------------------- if keyword_set(identi) then begin tmpstr = strmid(datset[0].datdir, 0, strlen(datset[0].datdir)-1) tmpstr = strjoin(strsplit(tmpstr, '-', /extract)) datset[i].identi = strjoin(strsplit(tmpstr, '-', /extract)) + $ '-' + string(byte(65+randomu(seed)*25)) + $ string(j, format='(I03)') endif ; INCREASE THE COUNTER FOR ELEMENTS IN ONE FOLDER ;----------------------------------------------------------------------- j += 1 endfor ; DESTROY THE CALIBRATOR DATABASE ;------------------------------------------------------------------------------- if keyword_set(caldat) then obj_destroy, boecal ; FILL IN DEFAULT VALUES FOR THE REDUCTION PARAMETERS ;------------------------------------------------------------------------------- if keyword_set(redpar) then begin datset.redpar.tsmoot = 10. datset.redpar.gsmoot = 20. datset.redpar.davera = 1 datset.redpar.skypos = [9.,25.] datset.redpar.nophot = 1 endif ; GUESS THE ASSOCIATED CALIBRATOR ;------------------------------------------------------------------------------- if keyword_set(caldat) then for i=0,n_elements(datset)-1 do begin penlty = abs(filtbl[i].mjdobs - filtbl.mjdobs) nocali = where(~ strmatch(datset.dprcat, '*CALIB*', /fold)) if nocali[0] ge 0 then nocali = [i, nocali] else nocali = i penlty[nocali] = !values.d_infinity tmpstr = min(penlty, calidx) datset[i].redpar.calibr = datset[calidx].identi endfor END ;+ ; NAME: ; MIDI_GETDAT ; ; PURPOSE: ; Read raw and calibrated MIDI data from EWS into a struct. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; result = midi_getdat(datset[, /rawdat] [, /fixbli] [, /silent]) ; ; REQUIRED INPUTS: ; datset a MIDI data set or list of MIDI data sets ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORDS: ; rawdat if set, the raw, uncalibrated data will be read ; fixbli fix the baseline info by calculating if from the uv coords ; silent suppress the output of informational messages ; ; OUTPUTS: ; result MIDI result structure containing the reduced data ; ; DESCRIPTION AND EXAMPLE: ; This programme reads MIDI data that has been reduced by EWS into a MIDI ; result struct. By default, calibrated data is read. If the keyword ; RAWDAT is set then the raw data will be read. The data is read using ; the oirGetVis and oirGetData routines on the appropriate files. Note ; that the programme can only read one type of data into an array of ; struct, e.g. only PRISM or GRISM data. Structures with mixed data types ; are not (yet) possible. The quality flags are supposed to hold the ; following information encoded in a byte value: ; position value meaning ; 0 1 data is bad ; 1 2 bad correlated flux ; 2 4 bad total flux ; 3 8 bad visibility ; 4 16 spectral calibrator ; 5 32 photometric calibrator ; 6 64 unused ; 7 128 unused ; The programme is called by: ; result = midi_getdat(sciset) ; ; CALLED BY: ; midi_redset ; midi_adu2jy ; ; CALLING: ; midi_getkey ; ; MODIFICATION HISTORY: ; 2008-08-26 Written by Konrad R. Tristram ; 2008-09-17 Konrad R. Tristram: Create named structures for PRISM/GRISM. ; 2009-03-13 Konrad R. Tristram: Added multiplication with COSPHI. ; 2009-05-21 Konrad R. Tristram: Added keyword FIXBLI. ; 2010-09-22 Konrad R. Tristram: Changed tag LSTIME to MJDOBS and double. ; 2011-08-18 Konrad R. Tristram: Catch data sets without photometry. ; 2012-02-07 Konrad R. Tristram: Added handling of incomplete data sets. ; 2012-06-12 Konrad R. Tristram: Added tags RASCEN, DECLIN and QUALIT. ; 2012-06-14 Konrad R. Tristram: Prevent crash for "corrupted" files. ; 2012-08-29 Konrad R. Tristram: Added more details on QUALIT tag. ; 2012-11-17 Konrad R. Tristram: Corrected error in calculation of PANGLE. ; FUNCTION MIDI_GETDAT, DATSET, RAWDAT=RAWDAT, FIXBLI=FIXBLI, SILENT=SILENT ; INITIALISE OR SET DEFAULT VALUES FOR A FEW VARIABLES ;------------------------------------------------------------------------------- numfil = n_elements(datset) ; SET SUFFIXES DEPENDING ON IF WE HAVE CALIBRATED OR UNCALIBRATED DATA ;------------------------------------------------------------------------------- if keyword_set(rawdat) then begin suffix = ['.corr', '.photometry', '.redcal'] endif else begin suffix = ['.calcorr', '.calphot', '.calvis'] endelse ; CHECK WHICH FILES ACTUALLY EXIST ;------------------------------------------------------------------------------- filist = transpose([[datset.ewstag + suffix[0] + '.fits'], $ [datset.ewstag + suffix[1] + '.fits'], $ [datset.ewstag + suffix[2] + '.fits']]) exists = file_test(filist) if total(exists) eq 0 then begin message, 'No data found! Exiting.', /continue return, -1 endif ; LOAD FIRST ELEMENT TO CHECK WHAT KIND OF DATA THIS IS ;------------------------------------------------------------------------------- tmpdat = (where(exists eq 1))[0] if (tmpdat mod 3) eq 1 then wavele = oirgetwavelength(filist[tmpdat]) $ else tmpdat = oirGetVis(filist[tmpdat], wave=wavele) ;wavele = oirgetwavelength(datset[0].filnam[0]) wavelm = n_elements(wavele) ; PRINT WHAT KIND OF DATA WILL BE LOADED ;------------------------------------------------------------------------------- if ~ keyword_set(silent) then begin case wavelm of 171: begin message, 'Reading PRISM data.', /continue strnam = 'midi_prism' end 151: begin message, 'Reading old PRISM data.', /continue strnam = 'midi_prism_old' end 261: begin message, 'Reading GRISM data.', /continue strnam = 'midi_grism' end else: message, 'Trying to read unknown data type.', /continue endcase endif ;CREATE A STRUCT HOLDING THE RESULT ;------------------------------------------------------------------------------- result = {identi:'', object:'', $ rascen:0.D, declin:0.D, $ pbline:0.0, pangle:0.0, $ ucoord:0.0, vcoord:0.0, $ baseli:'', mjdobs:0.D, $ wavele:fltarr(wavelm), $ 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), $ qualit:0B, commen:''} if n_elements(strnam) eq 1 then result = create_struct(name=strnam, result) result = replicate(result, numfil) ;COPY SOME INFORMATION FROM THE DATA SET STRUCT TO THE RESULT STRUCT ;------------------------------------------------------------------------------- result.identi = datset.identi result.object = datset.object result.commen = datset.commen ; SET THE HEADER KEYWORDS THAT NEED TO BE READ ;------------------------------------------------------------------------------- keywords = ['RA', 'DEC', 'MJD-OBS', 'ISS CONF STATION1', 'ISS CONF STATION2', $ 'ISS PBL12 END', 'ISS PBLA12 END', 'OBS TARG NAME'] ; LOOP OVER THE INDIVIDUAL DATASET ELEMENTS ;------------------------------------------------------------------------------- for i=0,numfil-1 do begin ; READ IN THE HEADER KEYWORDS FROM THE FIRST EXISTING FILE ;------------------------------------------------------------------------------- tmpdat = (where(exists[*,i] eq 1))[0] if tmpdat ge 0 then begin keytbl = midi_getkey(filist[tmpdat,i], keywords) ;COPY THE HEADER INFORMATION INTO THE STRUCT ;----------------------------------------------------------------------- result[i].rascen = keytbl.ra result[i].declin = keytbl.dec result[i].mjdobs = keytbl.mjdobs result[i].pbline = keytbl.isspbl12end result[i].pangle = keytbl.isspbla12end ; Note: ESO header keywords 'ISS PBL12 END' and 'ISS PBLA12 END' may ; be incorrect! Therefore the PA and BL should be calculated at the ; end from the UV coordinates determined by EWS. result[i].baseli = strtrim(keytbl.issconfstation1, 2) + '-' + $ strtrim(keytbl.issconfstation2, 2) ; Use the object name from the FITS file instead of from DATSET: ;result[i].object = keytbl.obstargname endif else result[i].qualit = 1B ; TRY TO LOAD THE INTERFEROMETRIC DATA FROM THE FILE ;------------------------------------------------------------------------------- if file_test(datset[i].ewstag + suffix[0] + '.fits') then begin tmpdat = oirGetVis(datset[i].ewstag + suffix[0] + '.fits', $ wave=wavele, ierr=ierror) endif else ierror = 1B ; CHECK IF INTERFEROMETRIC DATA WAS LOADED SUCCESSFULLY ;------------------------------------------------------------------------------- if ~ ierror then begin ; COPY THE UV COORDINATES INTO THE RESULT STRUCT ;----------------------------------------------------------------------- result[i].ucoord = tmpdat.ucoord result[i].vcoord = tmpdat.vcoord ; CALCULATE THE REDUCTION FACTOR FOR THE CORRELATED FLUX ;----------------------------------------------------------------------- cosphi = cos(tmpdat.visphi / 180. * !dpi) cosphi = 1.0 ; COPY THE CORRELATED FLUX AND PHASE DATA TO THE RESULT STRUCT ;----------------------------------------------------------------------- result[i].wavele = wavele result[i].corflx = tmpdat.visamp * cosphi result[i].corflxerr = tmpdat.visamperr * cosphi result[i].visphi = tmpdat.visphi result[i].visphierr = tmpdat.visphierr endif else result[i].qualit += 2B ; TRY TO LOAD THE PHOTOMETRY DATA FROM THE FILE ;------------------------------------------------------------------------------- if file_test(datset[i].ewstag + suffix[1] + '.fits') then begin tmpdat = oirGetData(datset[i].ewstag + suffix[1] + '.fits', ierr=ierror) endif else ierror = 1B ; CHECK IF PHOTOMETRY DATA WAS LOADED SUCCESSFULLY ;------------------------------------------------------------------------------- if ~ ierror then begin ; PHOTOMETRY HAS TO BE LOADED FROM DIFFERENT POSITIONS IN THE STRUCTS ;----------------------------------------------------------------------- if keyword_set(rawdat) then begin result[i].mskflx = tmpdat[5].data1 result[i].mskflxerr = tmpdat[11].data1 result[i].totflx = 0.5*(tmpdat[0].data1+tmpdat[1].data1) result[i].totflxerr = sqrt(tmpdat[6].data1^2+tmpdat[7].data1^2) endif else begin result[i].mskflx = tmpdat[1].data1 result[i].mskflxerr = tmpdat[3].data1 result[i].totflx = tmpdat[0].data1 result[i].totflxerr = tmpdat[2].data1 endelse ; IF WAVELENGTH STILL UNDEFINED THEN COPY IT NOW ;----------------------------------------------------------------------- if (total(result[i].wavele) eq 0) then begin wavele = datset[i].ewstag + suffix[1] + '.fits' result[i].wavele = oirGetWavelength(wavele) endif endif else result[i].qualit += 4B ; TRY TO LOAD THE VISIBILITY DATA FROM THE FILE ;------------------------------------------------------------------------------- if file_test(datset[i].ewstag + suffix[2] + '.fits') then begin tmpdat = oirGetVis(datset[i].ewstag + suffix[2] + '.fits', $ wave=wavele, ierr=ierror) endif else ierror = 1B ; CHECK IF VISIBILITY DATA WAS LOADED SUCCESSFULLY ;------------------------------------------------------------------------------- if ~ ierror then begin ; COPY THE UV COORDINATES INTO THE RESULT STRUCT IF NOT ALREADY SET ;----------------------------------------------------------------------- if (result[i].ucoord eq 0) and (result[i].vcoord eq 0) then begin result[i].ucoord = tmpdat.ucoord result[i].vcoord = tmpdat.vcoord endif ; COPY THE VISIBILITY DATA TO THE RESULT STRUCT ;----------------------------------------------------------------------- result[i].visamp = tmpdat.visamp * cosphi result[i].visamperr = tmpdat.visamperr * cosphi ; IF WAVELENGTH STILL UNDEFINED THEN COPY IT NOW ;----------------------------------------------------------------------- if (total(result[i].wavele) eq 0) then result[i].wavele = wavele endif else result[i].qualit += 8B ; END OF LOOP OVER THE INDIVIDUAL DATASET ELEMENTS ;------------------------------------------------------------------------------- endfor ; CALCULATE THE PROJECTED BASELINE LENGTH AND ANGLE FROM THE UV COORDINATES ;------------------------------------------------------------------------------- if keyword_set(fixbli) then begin result.pbline = sqrt(result.ucoord^2+result.vcoord^2) result.pangle = atan(result.ucoord, result.vcoord) / !dpi*180. endif ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, result END ;+ ; NAME: ; MIDI_CALINF ; ; PURPOSE: ; Get calibrator information for a list of calibrated MIDI results. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; calinf = midi_calinf(caldat, datset) ; ; REQUIRED INPUTS: ; caldat MIDI result structure containing calibrated data ; datset the MIDI data set used to reduce the data in CALDAT ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORDS: ; None. ; ; OUTPUTS: ; calinf structure containing information on the calibrators ; ; DESCRIPTION AND EXAMPLE: ; This programme looks for a data set in DATSET with an identity tag ; matching the identity tag in the result struct CALDAT. If a data set ; is found, then the calibrator for the data set is found in the same ; way, that is by searching the matching calibrator identity tag in the ; data set DATSET. If the corresponding calibrator was found, its flux ; and name is copied to the result struct. The result struct hence ; contains the index of the calibrated result and the index of the ; calibrator in DATSET as well as the flux and name of the calibrator. ; The program is called by: ; calinf = midi_calinf(caldat, datset) ; ; CALLED BY: ; None. ; ; CALLING: ; None. ; ; MODIFICATION HISTORY: ; 2008-08-26 Written by Konrad R. Tristram ; FUNCTION MIDI_CALINF, CALDAT, DATSET ; CREATE THE STRUCTURE HOLDING THE RESULT ;------------------------------------------------------------------------------- calinf = {sciidx:-1L, calidx:-1L, calnam:'', calflx:0.0} calinf = replicate(calinf, n_elements(caldat)) ; LOOP OVER ALL ELEMENTS IN THE CALIBRATED DATA RESULT ;------------------------------------------------------------------------------- for i=0,n_elements(caldat)-1 do begin $ ; FIND THE CORRESPONDING INDEX OF THE SCIENCE DATA IN THE DATSET ;---------------------------------------------------------------------- sciidx = where(caldat[i].identi eq datset.identi) calinf[i].sciidx = sciidx ; CHECK IF A VALID INDEX FOR THE SCIENCE DATA WAS FOUND ;---------------------------------------------------------------------- if sciidx ge 0 then begin ; FIND THE CORRESPONDING INDEX OF THE CALIBRATOR IN THE DATSET ;-------------------------------------------------------------- calidx = where(datset[sciidx].redpar.calibr eq datset.identi) calinf[i].calidx = calidx endif else begin ; PRINT AN ERROR MESSAGE AND SET CALIBRATOR INDEX TO -1 ;-------------------------------------------------------------- message, 'Could not find data set corresponding to res' + $ 'ult' + string(i, format='(I3)') + '.', /continue calinf[i].calidx = -1 endelse ; COPY CALIBRATOR NAME AND FLUX IF THE INDEX IS VALID ;---------------------------------------------------------------------- if calidx ge 0 then begin calinf[i].calnam = datset[calidx].object calinf[i].calflx = datset[calidx].calflx endif else begin message, 'Could not find calibrator corresponding to res' + $ 'ult' + string(i, format='(I3)') + '.', /continue endelse endfor ;RETURN THE RESULT ;------------------------------------------------------------------------------- return, calinf END ;+ ; NAME: ; MIDI_RVBCAL ; ; PURPOSE: ; Recalibrate fluxes with templates by Roy van Boekel - OBSOLETE ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; result = midi_rvbcal(caldat, calnam [, calflx] [, /output]) ; ; REQUIRED INPUTS: ; caldat MIDI result structure containing calibrated data ; calnam name of the calibrators used for the primary calibration ; ; OPTIONAL INPUTS: ; calflx flux of the calibrator that was used in midicalibrate ; ; KEYWORDS: ; output make a diagnostioc plot ; ; OUTPUTS: ; result MIDI result structure with recalibrated fluxes ; ; DESCRIPTION AND EXAMPLE: ; This function looks for a calibrators in the calibrator database by ; Roy van Boekel (which must be present in the !CALCAT system variable) ; matching the names in CALNAM. If no matching calibrator is found, the ; correction is skipped for the respective science source. If a ; calibrator is found, then a correction factor is calculated: it is the ; the fraction between the template spectrum and the approximation by a ; Rayleigh-Jeans law with a 10 micron flux as specified in CALFLX. If ; CALFLX is not specified, a flux value of 1 Jy will be assumed. This ; correction factor is then applied to the fluxes of the MIDI data. ; If no corresponding calibrator is found, the data is left unchanged. ; Optionally a plot showing the correction factor as well as the fluxes ; before and after the calibration can be generated. The program is ; called by: ; result = midi_rvbcal(caldat, calnam, calflx) ; ; CALLED BY: ; None. ; ; CALLING: ; None. ; ; MODIFICATION HISTORY: ; 2008-08-26 Written by Konrad R. Tristram ; 2009-05-12 Konrad R. Tristram: Modified to cope with new CALCAT struct. ; 2010-07-13 Konrad R. Tristram: Also use spline interpolation as in EWS. ; FUNCTION MIDI_RVBCAL, CALDAT, CALNAM, CALFLX, OUTPUT=OUTPUT ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- numele = n_elements(caldat) if n_elements(calnam) ne numele then begin message, 'Number of calibrator names does not match number of ' + $ 'data sets! Returning.' return, -1 endif if n_elements(calflx) ne numele then calflx = replicate(1.0, numele) ; CHECK IF SYSTEM VARIABLE CONTAINING THE CALIBRATOR DATABASE IS PRESENT ;------------------------------------------------------------------------------- defsysv, '!calcat', exists = calpre if calpre eq 0 then begin message, 'Calibrator database not found! The database has to be' + $ 'loaded into the system variable !CALCAT.', /continue return, -1 endif ; CREATE THE RESULT STRUCT BY COPYING THE INPUT DATA ;------------------------------------------------------------------------------- result = caldat ; LOOP OVER ALL THE ELEMENTS TO BE CALIBRATED ;------------------------------------------------------------------------------- for i=0,numele-1 do begin ; EXTRACT ONLY THE NUMBERS FROM THE CALIBRATOR NAME ;----------------------------------------------------------------------- calnum = byte(calnam[i]) numidx = where((calnum ge 48b) and (calnum le 57b)) calnum = strtrim(long(string(calnum[numidx])),2) ; FIND THE CALIBRATOR IN THE DATA BASE ;----------------------------------------------------------------------- calidx = where(!calcat.name eq 'hd_' + calnum) if calidx eq -1 then calidx = where(!calcat.name eq 'HD' + calnum) ; IF NOT FOUND THEN PRINT ERROR MESSAGE AND SKIP THE CORRECTION ;----------------------------------------------------------------------- if calidx eq -1 then begin message, 'The calibrator ' + calnam[i] + ' was not found. ' + $ 'Returning.', /continue ; IF THE CALIBRATOR WAS FOUND THEN APPLY THE CORRECTION ;----------------------------------------------------------------------- endif else begin print, i, calnam[i], format='("Data set ", I3, ": The ' + $ 'calibrator ", A, " was found.")' ; COMPARE FLUX GIVEN IN PROGRAMME CALL AND FLUX GIVEN IN THE DATABASE ; (THE DATABASE FLUX IS OBVIOUSLY FOR 10 MICRON AND NOT FOR 12 MICRON) ;----------------------------------------------------------------------- print, ' Flux for Rayleigh Jeans law is ' + $ string(calflx[i], format='(F6.3)') + ' Jy.' print, ' Flux from calibrator the database is ' + $ string(!calcat[calidx].f12.fnu, format='(F6.3)') + ' Jy.' print, ' Flux from spectrum in the database is ' + $ string(interpol(!calcat[calidx].spec.fnu, $ !calcat[calidx].spec.lam, 10), $ format='(F6.3)') + ' Jy.' ; CALCULATE THE CALIBRATOR FLUX FOR THE RESPECTIVE WAVELENGTH BINS ;----------------------------------------------------------------------- wavele = caldat[i].wavele caltem = interpol(!calcat[calidx].spec.fnu, !calcat[calidx].spec.lam, $ caldat[i].wavele, /spline) ; CALCULATE CORRECTION FACTOR BETWEEN RAYLEIGH JEANS AND THE TEMPLATE ;----------------------------------------------------------------------- correc = caltem / (calflx[i] * 100 / wavele^2) ; APPLY CORRECTION FACTOR TO THE PHOTOMETRIES AND THE CORRELATED FLUX ;----------------------------------------------------------------------- result[i].corflx = correc * result[i].corflx result[i].mskflx = correc * result[i].mskflx result[i].totflx = correc * result[i].totflx result[i].corflxerr = correc * result[i].corflxerr result[i].mskflxerr = correc * result[i].mskflxerr result[i].totflxerr = correc * result[i].totflxerr ; START OUTPUT OF DIAGNOSTIC PLOTS IF DESIRED ;----------------------------------------------------------------------- if keyword_set(output) then begin ; SAVE THE OLD PLOT SETTINGS, THEN CHANGE THE CHARACTER SIZE ;--------------------------------------------------------------- pltset = !p !p.charsize = 1.25 ; OPEN NEW PLOT WINDOW AND SET THE COLOURS ;--------------------------------------------------------------- window, /free, xsize=600, ysize=600 device, decompose=0 color1=165 color2=165 color3=200 color4=250 ; SELECT WAVELENGTH RANGE ;--------------------------------------------------------------- waveel = where((wavele gt 7.7) and (wavele lt 13.3)) ; DETERMINE SCALES FOR THE CORRECTION FACTOR PLOT ;--------------------------------------------------------------- yrange = 1.1 * max(abs(correc[waveel] - 1.)) yrange = [1-yrange, 1+yrange] ; MAKE CORRECTION FACTOR PLOT ;--------------------------------------------------------------- loadct, 39, /silent plot, [0,100],[1,1], xrange=[7.5, 13.5], yrange=yrange, $ xstyle=1, ystyle=0, xtitle='wavelength in !4l!Xm', $ ytitle='correction factor', xticklen=0.04, $ position=[60,50,570,220], /DEVICE loadct, 33, /silent oplot, wavele, correc, color=color1, thick=2 ; DETERMINE SCALES FOR THE FLUX PLOT ;--------------------------------------------------------------- flxdat = [result[i].corflx[waveel], $ result[i].mskflx[waveel], $ result[i].totflx[waveel], $ result[i].corflx[waveel] / correc, $ result[i].mskflx[waveel] / correc, $ result[i].totflx[waveel] / correc] yrange = 1.1 *max(flxdat) * [0.0, 1.1] ; MAKE FLUX PLOT ;--------------------------------------------------------------- loadct, 39, /silent plot, [0,100],[0,0], xrange=[7.5, 13.5], yrange=yrange, $ xstyle=1, ystyle=0, ytitle='flux in Jy', xticklen=0.02, $ xtickname=replicate(' ', 10), positi=[60,220,570,560], $ title='Recalibration with van Boekel template', $ /DEVICE, /NOERASE loadct, 33, /silent oplot, wavele, result[i].totflx, color=color2 oplot, wavele, caldat[i].totflx, color=color2, line=2 oplot, wavele, result[i].mskflx, color=color3 oplot, wavele, caldat[i].mskflx, color=color3, line=2 oplot, wavele, result[i].corflx, color=color4 oplot, wavele, caldat[i].corflx, color=color4, line=2 ; DETERMINE Y-RANGE FINALLY USED IN THE FLUX PLOT ;--------------------------------------------------------------- yrange = !y.crange[1] ; PLOT AND WRITE LEGEND ;--------------------------------------------------------------- oplot, [7.7, 8.0], 0.93 * yrange * [1,1], color=color2 oplot, [7.7, 8.0], 0.87 * yrange * [1,1], color=color2, line=2 oplot, [7.7, 8.0], 0.81 * yrange * [1,1], color=color3 oplot, [7.7, 8.0], 0.75 * yrange * [1,1], color=color3, line=2 oplot, [7.7, 8.0], 0.69 * yrange * [1,1], color=color4 oplot, [7.7, 8.0], 0.63 * yrange * [1,1], color=color4, line=2 loadct, 39, /silent xyouts, 8.1, 0.92 * yrange, 'new total flux' xyouts, 8.1, 0.86 * yrange, 'old total flux' xyouts, 8.1, 0.80 * yrange, 'new masked flux' xyouts, 8.1, 0.74 * yrange, 'old masked flux' xyouts, 8.1, 0.68 * yrange, 'new correlated flux' xyouts, 8.1, 0.62 * yrange, 'old correlated flux' ; RESTORE THE ORIGINAL PLOT SETTINGS ;--------------------------------------------------------------- !p = pltset endif ; END OF THE CORRECTION PART (ONLY EXECUTED IF CALIBRATOR IS FOUND) ;----------------------------------------------------------------------- endelse endfor ;RETURN THE RESULT ;------------------------------------------------------------------------------- return, result END ;+ ; NAME: ; MIDI_REDPLT ; ; PURPOSE: ; Plot the fluxes and visibilities of a MIDI result structure ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; midi_redplt, result [, ctable] [, colour] [, xrange=[7.5, 13.5]] ; [, /rawdat] [, /plterr] [, psfile=psfile] ; ; REQUIRED INPUTS: ; result MIDI result structure ; ; OPTIONAL INPUTS: ; ctable the colour table to be used ; colour the colour indices for the different data sets ; ; KEYWORDS: ; xrange a two element array specifying the x-range ; rawdat indicates that raw data is plotted ; plterr also plot error bars for the data plotted ; psfile selects output to PS and specifies the file name ; ; OUTPUTS: ; None. ; ; DESCRIPTION AND EXAMPLE: ; This programmes plots the masked, total and correlated fluxes as well ; as the visibility and the phase into one single plot. Additionally the ; UV position of the observation is also plotted. By default axis labels ; are for calibrated data. For uncalibrated data, the corect axis labels ; are obtained by setting the keyword RAWDAT. The program is called by: ; midi_redplt, result ; ; CALLED BY: ; midi_redset ; ; CALLING: ; None. ; ; MODIFICATION HISTORY: ; 2008-08-26 Written by Konrad R. Tristram ; 2009-05-12 Konrad R. Tristram: Added possibility to plot error bars. ; 2010-08-07 Konrad R. Tristram: Changed the way the xrange is set. ; 2011-08-12 Konrad R. Tristram: Added keyword XRANGE. ; 2012-10-15 Konrad R. Tristram: Change thickness when plotting errors. ; 2012-11-01 Konrad R. Tristram: Adjust PS font size for !p.font=0. ; PRO MIDI_REDPLT, RESULT, CTABLE, COLOUR, XRANGE=XRANGE, RAWDAT=RAWDAT, $ PLTERR=PLTERR, PSFILE=PSFILE ; INITIALISE OR SET DEFAULT VALUES FOR A FEW VARIALBLES ;------------------------------------------------------------------------------- if n_elements(xrange) ne 2 then xrange = [7.5, 13.5] erroff = 0.01 ; offset of error bars errwid = 0.00 ; width of error bar ends numdat = n_elements(result) if keyword_set(plterr) then plterr = 1 else plterr = 0 if n_elements(colour) ne numdat then begin colour = 165 + (indgen(numdat) / ((numdat - 1.) > 1.))^1.5 * (254-165) ctable = 33 endif ; PREPARE VARIABLES DEPENDING ON IF THIS IS CALIBRATED OR UNCALIBRATED DATA ;------------------------------------------------------------------------------- if keyword_set(rawdat) then begin plttit = ['Raw masked flux (masked photometry, geometric mean) ', $ 'Raw total flux (total photometry, arithmetric mean) ', $ 'Raw correlated flux', 'Raw visibility', $ 'Raw differential phase'] ytitle = ['raw masked flux F!Dmsk,raw!N [ADU]', $ 'raw total flux F!Dtot,raw!N [ADU]', $ 'raw correlated flux F!Dcor,raw!N [ADU]', $ 'raw visibility V', 'raw differential phase' + $ '!4u!X!Ddif,raw!N [' + string("260B) + ']'] endif else begin plttit = ['Calibrated masked flux (masked phot., geometric mean) ', $ 'Calibrated total flux (total phot., arithmetric mean) ', $ 'Calibrated correlated flux', 'Calibrated visibility', $ 'Calibrated differential phase'] ytitle = ['masked flux F!Dmsk!N [Jy]', 'total flux F!Dtot!N [Jy]', $ 'correlated flux F!Dcor!N [Jy]', 'visibility V', $ 'differential phase !4u!X!Ddif!N [' + string("260B) + ']'] endelse ; MAKE A BACKUP OF THE PLOT SETTINGS ;------------------------------------------------------------------------------- p=!p !p.multi=[0,2,3] d_old=!d.name ; OPEN A NEW PLOT WINDOW OR FILE AND SET THE PLOT ENVIRONMENT ;------------------------------------------------------------------------------- if keyword_set(psfile) then begin set_plot, 'PS' device, file=psfile + '.ps', /portrait, color=1, $ xoffset=1.0, yoffset=3.3, xsize=19., ysize=24. if !p.font eq 0 then !p.charsize=0.65 else !p.charsize=0.75 position = [14500, 600, 18000, 7600] legende = [9800, 600, 400] endif else begin set_plot, 'X' device, decompose=0 window, xsize=800, ysize=1000, /FREE !p.charsize=1.25 position = [600, 020, 750, 320] legende = [420, 30, 15] endelse ; SELECT WAVELENGTH-ELEMENTS BETWEEN 8 AND 13 MICRON (FOR SCALING THE GRAPHS) ;------------------------------------------------------------------------------- wavsel = where(result.wavele gt 8.0 and result.wavele lt 13.0) ; PLOT THE MASKED PHOTOMETRY ;------------------------------------------------------------------------------- loadct, 0, /silent ;minplt = min((result.mskflx)[wavsel]) ;maxplt = max((result.mskflx)[wavsel]) minplt = min([(result.mskflx)[wavsel], (result.totflx)[wavsel]]) maxplt = max([(result.mskflx)[wavsel], (result.totflx)[wavsel]]) ;if not keyword_set(rawdat) then maxplt = maxplt < 40 plot, xrange, [minplt, maxplt], xstyle=1, ystyle=0, title=plttit[0], $ xmargin=[10, 2], xtitle='wavelength !4k!X [!4l!Xm]', ytitle=ytitle[0], $ charsize=1.75*!p.charsize, /NODATA if !y.crange[0] lt 0.0 then oplot, [0,100], [0,0], color=175 loadct, ctable, /silent for i=0,numdat-1 do begin oplot, result[i].wavele, result[i].mskflx, color=colour[i], $ thick=plterr+1 if plterr then errplot, result[i].wavele+erroff*i, result[i].mskflx-$ result[i].mskflxerr, result[i].mskflx+$ result[i].mskflxerr, color=colour[i], width=errwid endfor ; PLOT THE UNMASKED PHOTOMETRY ;------------------------------------------------------------------------------- loadct, 0, /silent ;minplt = min((result.totflx)[wavsel]) ;maxplt = max((result.totflx)[wavsel]) plot, xrange, [minplt, maxplt], xstyle=1, ystyle=0, title=plttit[1], $ xmargin=[10, 2], xtitle='wavelength !4k!X [!4l!Xm]', ytitle=ytitle[1], $ charsize=1.75*!p.charsize, /NODATA if !y.crange[0] lt 0.0 then oplot, [0,100], [0,0], color=175 loadct, ctable, /silent for i=0,numdat-1 do begin oplot, result[i].wavele, result[i].totflx, color=colour[i], $ thick=plterr+1 if plterr then errplot, result[i].wavele+erroff*i, result[i].totflx-$ result[i].totflxerr, result[i].totflx+$ result[i].totflxerr, color=colour[i], width=errwid endfor ; PLOT THE CORRELATED FLUX ;------------------------------------------------------------------------------- loadct, 0, /silent minplt = 0 maxplt = max((result.corflx)[wavsel]) ;if not keyword_set(rawdat) then maxplt = 3.0 plot, xrange, [minplt, maxplt], xstyle=1, ystyle=0, title=plttit[2], $ xmargin=[10, 2], xtitle='wavelength !4k!X [!4l!Xm]', ytitle=ytitle[2], $ charsize=1.75*!p.charsize, /NODATA loadct, ctable, /silent for i=0,numdat-1 do begin oplot, result[i].wavele, result[i].corflx, color=colour[i], $ thick=plterr+1 if plterr then errplot, result[i].wavele+erroff*i, result[i].corflx-$ result[i].corflxerr, result[i].corflx+$ result[i].corflxerr, color=colour[i], width=errwid endfor ; PLOT THE VISIBILITY ;------------------------------------------------------------------------------- loadct, 0, /silent minplt = 0 maxplt = max((result.visamp)[wavsel]) < 5.0 plot, xrange, [minplt, maxplt], xstyle=1, ystyle=0, title=plttit[3], $ xmargin=[10, 2], xtitle='wavelength !4k!X [!4l!Xm]', ytitle=ytitle[3], $ charsize=1.75*!p.charsize, /NODATA if !y.crange[1] gt 1.0 then oplot, [0,100], [1,1], color=175 loadct, ctable, /silent for i=0,numdat-1 do begin oplot, result[i].wavele, result[i].visamp, color=colour[i], $ thick=plterr+1 if plterr then errplot, result[i].wavele+erroff*i, result[i].visamp-$ result[i].visamperr, result[i].visamp+$ result[i].visamperr, color=colour[i], width=errwid endfor ; PLOT THE DIFFERENTIAL PHASE ;------------------------------------------------------------------------------- loadct, 0, /silent minplt = min((result.visphi)[wavsel]) > (-180) maxplt = max((result.visphi)[wavsel]) < (+180) plot, xrange, [minplt, maxplt], xstyle=1, ystyle=0, title=plttit[4], $ xmargin=[10, 2], xtitle='wavelength !4k!X [!4l!Xm]', ytitle=ytitle[4], $ charsize=1.75*!p.charsize, /NODATA oplot, [0,100], [0,0], color=175 loadct, ctable, /silent for i=0,numdat-1 do begin oplot, result[i].wavele, result[i].visphi, color=colour[i], $ thick=plterr+1 if plterr then errplot, result[i].wavele+erroff*i, result[i].visphi-$ result[i].visphierr, result[i].visphi+$ result[i].visphierr, color=colour[i], width=errwid endfor ; PLOT A SMALL UV PLANE WITH THE LOCATION OF THE MEASUREMENTS ;------------------------------------------------------------------------------- uvsize = 1.05*max(sqrt(result.ucoord^2 + result.vcoord^2)) loadct, 0, /silent plot, [uvsize, 0.0001], [-uvsize, uvsize], xstyle=5, ystyle=5, $ xrange=[uvsize, 0.0001], xmargin=[0,0], ymargin=[0,0], $ position=position, /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 axis, 0 , 0, /xaxis, xstyle=1, xtickinterval=fix(uvsize/10-0.5)*5, $ xtickname=replicate(' ', 10), ticklen=0.02, charsize=1.75*!p.charsize axis, 0 , 0, /yaxis, ystyle=1, ytickinterval=fix(uvsize/10-0.5)*5, $ ticklen=0.04, charsize=1.75*!p.charsize xyouts, 0.07 * uvsize, 1.03 * uvsize, 'V' xyouts, 1.07 * uvsize, -0.07 * uvsize, 'U' loadct, ctable, /silent for i=0,numdat-1 do begin if result[i].ucoord lt 0.0 then begin xyouts, -result[i].ucoord, -result[i].vcoord, $ string(i, format='(I02)'), color=colour[i] endif else begin xyouts, result[i].ucoord, result[i].vcoord, $ string(i, format='(I02)'), color=colour[i] endelse endfor ; PLOT THE LEGEND ;------------------------------------------------------------------------------- loadct, ctable, /silent for i=0,numdat-1 do xyouts, legende[0], legende[1] + (numdat-i) * legende[2], $ string(i, format='(I02)') + ': ' + $ strmid(date_conv(result[i].mjdobs+2400000.5, 'F'), 0, 16) + ' ' + $ ; '20'+strjoin(strmid(result[i].identi, [0,2,4], 2), '-') + ' ' + $ result[i].object, color=colour[i], $ charsize=1.1*!p.charsize, /DEVICE ; RESET PLOT SETTINGS ;------------------------------------------------------------------------------- if keyword_set(psfile) then begin device, /close endif loadct, 39, /silent !p = p set_plot, d_old END ;+ ; NAME: ; MIDI_MODPHT ; ; PURPOSE: ; Modify the sky bands and recompile the photometry programmes of EWS. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; midi_modpht, skypos, (/verbose) ; ; REQUIRED INPUTS: ; skypos position of the sky bands ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; verbose print out the changes which were applied ; ; OUTPUTS: ; oirChopPhotoImagesMod.c and a modified oirChopPhotoImages executable ; ; DESCRIPTION AND EXAMPLE: ; This programme copies the source code from oirChopPhotoImages.c ; and oirMakePhotoSpectra.c into oirChopPhotoImagesMod.c and ; oirMakePhotoSpectraMod.c respectively in order to then change the ; definitions of the sky band positions. This is achieved by looping ; through the individual lines of the original file and altering the ; relevant lines before writing them to the output file. The resulting ; source code in the files named *Mod.c is compiled to create updated ; versions of the executables. The program is called by: ; midi_modpht, [6,22], /verbose ; ; CALLED BY: ; midi_redset ; ; CALLING: ; None. ; ; MODIFICATION HISTORY: ; 2010-07-12 Written by Konrad R. Tristram ; PRO MIDI_MODPHT, SKYPOS, VERBOSE=VERBOSE ; INITIALISE A FEW VARIABLES ;------------------------------------------------------------------------------- if keyword_set(verbose) then verbos = verbose else verbos = 0 strlin = '' ; CHANGE INTO THE DIRECTORY OF THE EWS PROGRAMMES ;------------------------------------------------------------------------------- cd, getenv('drsRoot') + '/c/src/', curren=curdir ; OPEN oirChopPhotoImages FOR READING AND AN OUTPUT FILE FOR WRITING ;------------------------------------------------------------------------------- openr, lun_in, 'exec/oirChopPhotoImages.c', /get_lun openw, lun_ou, 'exec/oirChopPhotoImagesMod.c', /get_lun ; LOOP OVER ALL LINES OF THE INPUT FILE ;------------------------------------------------------------------------------- while ~ eof(lun_in) do begin ; READ THE CURRENT INPUT LINE ;----------------------------------------------------------------------- readf, lun_in, strlin ; MODIFY THE STRING AND PRINT A MESSAGE FOR THE PRISMBOTTOM DEFINITION ;----------------------------------------------------------------------- if strmatch(strlin, '#define PRISMBOTTOM*') then begin tmpstr = string(skypos[0], format='(I3)') strlin = '#define PRISMBOTTOM' + tmpstr if verbos then print, 'Changed PRISMBOTTOM to' + tmpstr + '.' ; MODIFY THE STRING AND PRINT A MESSAGE FOR THE PRISMTOP DEFINITION ;----------------------------------------------------------------------- endif else if strmatch(strlin, '#define PRISMTOP*') then begin tmpstr = string(skypos[1], format='(I3)') strlin = '#define PRISMTOP' + tmpstr if verbos then print, 'Changed PRISMTOP to' + tmpstr + '.' endif ; WRITE THE CURRENT LINE INTO THE NEW FILE ;----------------------------------------------------------------------- printf, lun_ou, strlin endwhile ; CLOSE THE INPUT AND OUTPUT FILES ;------------------------------------------------------------------------------- close, lun_in, lun_ou ; RECOMPILE oirChopPhotoImages USING THE MODIFIED VERSION ;------------------------------------------------------------------------------- comman = 'gcc -o ../bin/oirChopPhotoImages -W -pedantic -I../../cfitsio ' + $ '-I../include -L../../cfitsio -L../lib -L../object ' + $ 'exec/oirChopPhotoImagesMod.c -lEWS -lcfitsio -lm' if not verbos then comman += ' -w' spawn, comman, exit_status=status ; OPEN oirMakePhotoSpectra FOR READING AND AN OUTPUT FILE FOR WRITING ;------------------------------------------------------------------------------- openr, lun_in, 'exec/oirMakePhotoSpectra.c' openw, lun_ou, 'exec/oirMakePhotoSpectraMod.c' ; LOOP OVER ALL LINES OF THE INPUT FILE ;------------------------------------------------------------------------------- while ~ eof(lun_in) do begin ; READ THE CURRENT INPUT LINE ;----------------------------------------------------------------------- readf, lun_in, strlin ; MODIFY THE STRING AND PRINT A MESSAGE FOR THE PRISMBOTTOM DEFINITION ;----------------------------------------------------------------------- if strmatch(strlin, ' int yBottom =*') then begin tmpstr = string(skypos[0], format='(I3)') strlin = ' int yBottom =' + tmpstr + ';' if verbos then print, 'Changed yBottom to' + tmpstr + '.' ; MODIFY THE STRING AND PRINT A MESSAGE FOR THE PRISMTOP DEFINITION ;----------------------------------------------------------------------- endif else if strmatch(strlin, ' int yTop =*') then begin tmpstr = string(skypos[1], format='(I3)') strlin = ' int yTop =' + tmpstr + ';' if verbos then print, 'Changed yTop to' + tmpstr + '.' endif ; WRITE THE CURRENT LINE INTO THE NEW FILE ;----------------------------------------------------------------------- printf, lun_ou, strlin endwhile ; CLOSE THE INPUT AND OUTPUT FILES ;------------------------------------------------------------------------------- close, lun_in, lun_ou free_lun, lun_in, lun_ou ; RECOMPILE oirMakePhotoSpectra USING THE MODIFIED VERSION ;------------------------------------------------------------------------------- comman = 'gcc -o ../bin/oirMakePhotoSpectra -W -pedantic -I../../cfitsio ' + $ '-I../include -L../../cfitsio -L../lib -L../object ' + $ 'exec/oirMakePhotoSpectraMod.c -lEWS -lcfitsio -lm' if not verbos then comman += ' -w' spawn, comman, exit_status=status ; RETURN TO THE CURRENT DIRECTORY ;------------------------------------------------------------------------------- cd, curdir END ;+ ; NAME: ; MIDI_REDSET ; ; PURPOSE: ; Fully reduce the data in a MIDI data set. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; midi_redset, datset [, prefix] [, /noredu] [, /noplot] [, /shimsk] ; [, /fringe] [, maxopd=maxopd] [, slopef=slopef] ; ; REQUIRED INPUTS: ; datset a MIDI data set or list of MIDI data sets ; ; OPTIONAL INPUTS: ; prefix a prefix for the PS files containing the result plots ; ; KEYWORDS: ; noredu skip the number crunching part of the reduction ; shimsk shift the mask ; fringe also calculate the fringe image from the data ; maxopd use a certain maxopd value ; slopef use the slope fitting routines ('T'-routines) ; ; OUTPUTS: ; three plots generated with MIDI_REDPLT: sciraw.ps, calraw.ps, scical.ps ; ; DESCRIPTION AND EXAMPLE: ; This programme loops through all elements in DATSET and reduces the ; data using 'midipipe' or 'faintpipe', or the respective '*vispipe' ; or '*photopipe' if the data set is incomplete. The pipeline routines ; are called using the parameters stored in the .REDPAR tag of DATSET. ; This first part of the reduction, which is the most time consuming, ; can be skipped by setting the keyword NOREDU. After the basic ; reduction, the data is calibrated with 'midicalibrate' using the ; calibrator database included in EWS. When no associated calibrator ; information is found, the calibration is skipped automatically. ; Finally three PS files are generated if the optional parameter ; PREFIX is specified. The files contain plots produced with ; MIDI_REDPLT of (1) the raw science data, (2) the raw calibrator ; data and (3) the calibrated science data, respectively. ; The programme is called by: ; midi_redset, datset ; ; CALLED BY: ; None. ; ; CALLING: ; midi_getdat, midi_redplt, midi_modpht ; ; MODIFICATION HISTORY: ; 2008-08-26 Written by Konrad R. Tristram ; 2009-05-20 Konrad R. Tristram: Added keyword SHIMSK for shifting masks. ; 2010-02-08 Konrad R. Tristram: Call midicalibrate twice to handle nophot ; 2010-02-08 Konrad R. Tristram: Added SKYMASK and changed mask handling ; 2010-04-14 Konrad R. Tristram: Added the PREFIX keyword ; 2010-07-12 Konrad R. Tristram: Added automatic adaptation of sky bands. ; 2010-07-13 Konrad R. Tristram: Added use of new calibration routines. ; 2010-08-12 Konrad R. Tristram: Added keyword FRINGE and MAXOPD ; 2010-11-18 Konrad R. Tristram: Added keyword SLOPEF ; 2011-08-18 Konrad R. Tristram: Adapted for reduction without photometry. ; 2012-02-07 Konrad R. Tristram: Improved handling of incomplete data. ; 2012-08-03 Konrad R. Tristram: Improved use of DETDIT for slope fitting. ; 2012-08-30 Konrad R. Tristram: Also use calibrators not in database. ; PRO MIDI_REDSET, DATSET, PREFIX, NOREDU=NOREDU, SHIMSK=SHIMSK, FRINGE=FRINGE, $ MAXOPD=MAXOPD, MSMOOTH=MSMOOTH, SLOPEF=SLOPEF ; IF SLOPEFITTING ROUTINES ARE NOT EXPLICITELY TURNED OFF, THEN USE THEM ;------------------------------------------------------------------------------- if n_elements(slopef) lt 1 then slopef = 1 if n_elements(msmooth) lt 1 then msmooth = 1.0 ; DETERMINE THE NUMBER OF ELEMENTS IN THE DATA SET ;------------------------------------------------------------------------------- numdat = n_elements(datset) - 1 ; LOAD THE RAINBOW COLOUR TABLE ;------------------------------------------------------------------------------- loadct, 39, /silent ; LOOP OVER ALL THE ELEMENTS IN THE DATA SET FOR THE MAIN DATA REDUCTION ;------------------------------------------------------------------------------- for i=0,numdat do begin ; TEST IF REDUCTION PATH EXISTS, OTHERWISE CREATE IT ;----------------------------------------------------------------------- ewsdir = path_sep() + 'ews' if file_test(datset[i].datdir + ewsdir, /directory) then begin ;print, 'EWS directory already exists. Using it.' endif else begin message, 'EWS directory not existing. Creating it.', /continue file_mkdir, datset[i].datdir + ewsdir endelse ; PRINT OUT A DIAGNOSTIC MESSAGE ;----------------------------------------------------------------------- message, string(i, datset[i].redpar.tsmoot, datset[i].redpar.gsmoot, $ datset[i].redpar.davera, datset[i].redpar.skypos[0], $ datset[i].redpar.skypos[1], $ format='("Data set", I3, ": Reducing data with ' + $ 'smo =", F5.1, ", gsm =", F6.2, ", dAve =", I2, ' + $ '", sky = [", I3, ",", I3, "].")'), /continue ; CHECK WHICH OF THE FILES ACTUALLY EXIST ;----------------------------------------------------------------------- exists = [(strsplit(datset[i].filnam[0], /extract))[0], $ (strsplit(datset[i].filnam[1], /extract))[0], $ (strsplit(datset[i].filnam[2], /extract))[0]] exists = file_test(exists) ; IF KEYWORD 'NORED' SET, PRETEND NO FILES EXIST TO SKIP THE REDUCTION ;----------------------------------------------------------------------- if keyword_set(noredu) then exists = [0,0,0] ; TEST WHETHER THE OBJECT MASK EXISTS ;----------------------------------------------------------------------- mskfil = datset[i].redpar.mskfil if ~ file_test(mskfil) then begin dummy = temporary(mskfil) message, 'Object mask not found or not provided. ' + $ 'Using standard EWS masks.', /continue endif ; TEST WHETHER THE SKY MASK EXITS ;----------------------------------------------------------------------- skyfil = datset[i].redpar.skyfil if ~ file_test(skyfil) then begin dummy = temporary(skyfil) message, 'Sky mask not found or not provided. ' + $ 'Using standard EWS sky bands.', /continue endif ; DETERMINE THE DIT AS A FACTOR FOR THE SMOOTHING PARAMETERS ;----------------------------------------------------------------------- ; Also assume that if gsmooth is less than 1.0, then it is in seconds if exists[0] and (slopef ge 1) and $ datset[i].redpar.gsmoot gt 1.0 then begin detdit = midiGetKeyword('DET DIT',datset[i].filnam[0]) endif else detdit = 1.0 ; IF ALL THREE FILES ARE PRESENT THEN RUN MIDIPIPE/FAINTPIPE ;----------------------------------------------------------------------- if total(exists) eq 3 then begin ; SET THE SKY BANDS AND RECOMPILE THE RELEVANT PROGRAMMES ;--------------------------------------------------------------- midi_modpht, datset[i].redpar.skypos ; RUN THE FAINT ROUTINE ;--------------------------------------------------------------- if slopef eq 1 then begin midipipe, datset[i].ewstag, datset[i].filnam, $ smooth=datset[i].redpar.tsmoot*detdit, $ gsmooth=datset[i].redpar.gsmoot*detdit, $ ave=datset[i].redpar.davera-1, $ mask=mskfil, minopd=minopd, $ maxopd=maxopd, msmooth=msmooth, twopass=1, $ dsky=datset[i].redpar.skydis, skyMask=skyfil ; RUN THE STANDARD ROUTINE ;--------------------------------------------------------------- endif else begin oldpipe, datset[i].ewstag, datset[i].filnam, $ smooth=datset[i].redpar.tsmoot*detdit, $ gsmooth=datset[i].redpar.gsmoot*detdit, $ dave=datset[i].redpar.davera, $ mask=mskfil, minopd=minopd, maxopd=maxopd, $ fast=slopef, shiftmask=shimsk, $ dsky=datset[i].redpar.skydis;, skyMask=skyfil endelse ; IF FRINGE TRACK IS PRESENT THEN RUN MIDI-/FAINTVISPIPE ;----------------------------------------------------------------------- endif else if exists[0] eq 1 then begin ; RUN THE FAINT ROUTINE ;--------------------------------------------------------------- if slopef then begin midiVisPipe, datset[i].ewstag, datset[i].filnam[0], $ smooth=datset[i].redpar.tsmoot*detdit, $ gsmooth=datset[i].redpar.gsmoot*detdit, $ ave=datset[i].redpar.davera-1, $ mask=mskfil, minopd=minopd, $ maxopd=maxopd, msmooth=msmooth, twopass=1 ; RUN THE STANDARD ROUTINE ;--------------------------------------------------------------- endif else begin oldVisPipe, datset[i].ewstag, datset[i].filnam[0], $ smooth=datset[i].redpar.tsmoot, $ gsmooth=datset[i].redpar.gsmoot*detdit, $ dave=datset[i].redpar.davera, $ mask=mskfil, minopd=minopd, $ maxopd=maxopd;, shiftmask=shimsk endelse ; IF BOTH PHOTOMETRIES ARE PRESENT THEN RUN MIDI-/FAINTPHOTOPIPE ;----------------------------------------------------------------------- endif else if array_equal(exists, [0,1,1]) then begin ; SET THE SKY BANDS AND RECOMPILE THE RELEVANT PROGRAMMES ;--------------------------------------------------------------- midi_modpht, datset[i].redpar.skypos ; RUN THE FAINT ROUTINE ;--------------------------------------------------------------- if slopef then begin midiPhotoPipe, datset[i].ewstag, mask=mskfil, $ datset[i].filnam[1:2], $ dsky=datset[i].redpar.skydis, $ skyMask=skyfil ; RUN THE STANDARD ROUTINE ;--------------------------------------------------------------- endif else begin midiPhotoPipe, datset[i].ewstag, mask=mskfil, $ datset[i].filnam[1:2], $ dsky=datset[i].redpar.skydis endelse endif else begin message, 'Skipping actual data reduction.', /continue endelse ; CHECK IF FRINGE TRACK IS PRESENT AND FRINGE IMAGE IS DESIRED ;----------------------------------------------------------------------- if (exists[0] eq 1) and keyword_set(fringe) then begin ; CHECK IF A VALID CORRELATED FLUX FILE WAS CREATED ;--------------------------------------------------------------- ; Only then it makes sense to also calculate the fringe image. detdit = oirgetvis(datset[i].ewstag+'.corr.fits', ierr=ierror) if ~ ierror then begin ; ONLY CALCULATE THE FRINGE IMAGE ;------------------------------------------------------- detdit = midiGetKeyword('DET DIT',datset[i].filnam[0]) midiFringeImage, datset[i].ewstag, datset[i].filnam[0],$ smooth=datset[i].redpar.tsmoot*detdit,$ noAve=datset[i].redpar.davera-1 endif endif endfor ; LOAD THE CALIBRATOR DATA BASE OF ROY VAN BOEKEL ;------------------------------------------------------------------------------- boecal = vboekelbase() ; LOOP AGAIN OVER ALL THE ELEMENTS IN THE DATA SET AND CALIBRATE THE DATA ;------------------------------------------------------------------------------- for i=0,numdat do begin ; FIND THE ASSOCIATED CALIBRATOR ;----------------------------------------------------------------------- asocal = (where(datset[i].redpar.calibr eq datset.identi))[0] ; CHECK IF AN ASSOCIATED CALIBRATOR WAS FOUND ;----------------------------------------------------------------------- if asocal ne -1 then begin ; PRINT OUT A DIAGNOSTIC MESSAGE ;--------------------------------------------------------------- caltag = strsplit(datset[asocal].ewstag, path_sep(), /extract) caltag = caltag[n_elements(caltag)-1] message, string(i, caltag, datset[asocal].calflx, $ datset[asocal].caldia, datset[i].redpar.nophot, $ format='("Data set", I3, ": Calibrating with ", ' + $ 'A, " (F =", F5.1, " Jy, d =", F5.1, ' + $ '" mas) and noPhot =", I2, ".")'), /continue ; CHECK IF THE NECESSARY FILES ACTUALLY EXIST ;--------------------------------------------------------------- docali = [file_test(datset[i].ewstag+'.corr.fits') and $ file_test(datset[asocal].ewstag+'.corr.fits'), $ file_test(datset[i].ewstag+'.photometry.fits') and $ file_test(datset[asocal].ewstag+'.photometry.fits')] ; CHECK IF CALIBRATOR IS IN THE DATA BASE ;--------------------------------------------------------------- srcdat = boecal->sourceFromFile(datset[asocal].ewstag + $ '.corr.fits', rad=0.02) ; IF CALIBRATOR NOT IN THE DATA BASE USE R.-J. APPROXIMATION ;--------------------------------------------------------------- if (size(srcdat, /type) ne 8) and $ finite(datset[asocal].calflx) and $ finite(datset[asocal].caldia) then calbas = 0B $ else calbas = boecal ; CALIBRATE THE SCIENCE DATA USING THE CALIBRATOR DATA ;--------------------------------------------------------------- if docali[0] and docali[1] then $ midicalibrate, datset[i].ewstag, datset[asocal].ewstag,$ calFlux10=datset[asocal].calflx, $ diam=datset[asocal].caldia, $ calDatabase=calbas if docali[0] and datset[i].redpar.nophot then $ midicalibrate, datset[i].ewstag, datset[asocal].ewstag,$ calFlux10=datset[asocal].calflx, $ diam=datset[asocal].caldia, $ calDatabase=calbas, /nophot endif endfor ;!x.style=0 ; This should not be necessary any more with newer than ;!y.style=0 ; August 2010 EWS versions. ;!x.range=0 ; DESTROY THE CALIBRATOR DATABASE ;------------------------------------------------------------------------------- obj_destroy, boecal ; START PLOTTING IF PLOTS ARE NOT TURNED OFF ;------------------------------------------------------------------------------- if keyword_set(prefix) then begin ; GET THE INDICES OF THE SCIENCE AND THE CALIBRATOR DATA ;----------------------------------------------------------------------- sciidx = where(datset.dprcat eq 'SCIENCE ') calidx = where(datset.dprcat eq 'CALIB ') ; PLOT THE RAW SCIENCE DATA ;----------------------------------------------------------------------- if sciidx[0] ge 0 then begin datres = midi_getdat(datset[sciidx], /rawdat, /silent) numdat = n_elements(datres) ctable = 33 colour = 164 + (indgen(numdat) / (numdat - 1.))^1.5 * (254-164) midi_redplt, datres, ctable, colour, /rawdat, $ psfile=prefix+'_sciraw' endif ; PLOT THE RAW CALIBRATOR DATA ;----------------------------------------------------------------------- if calidx[0] ge 0 then begin datres = midi_getdat(datset[calidx], /rawdat, /silent) numdat = n_elements(datres) ctable = 39 colour = 20 + (indgen(numdat) / (numdat - 1.))^1.0 * (120-20) midi_redplt, datres, ctable, colour, /rawdat, $ psfile=prefix+'_calraw' endif ; PLOT THE CALIBRATED SCIENCE DATA ;----------------------------------------------------------------------- if sciidx[0] ge 0 then begin datres = midi_getdat(datset[sciidx], /silent) numdat = n_elements(datres) ctable = 33 colour = 164 + (indgen(numdat) / (numdat - 1.))^1.5 * (254-164) midi_redplt, datres, ctable, colour, $ psfile=prefix+'_scical' endif endif END ;+ ; NAME: ; MIDI_PLTPRP ; ; PURPOSE: ; Plot a FITS header keyword for several MIDI data files. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; midi_pltprp, fileli, keywrd [, ytitle] [, psfile=psfile] ; [, yrange=yrange] [, titles=titles] [, labels=labels] ; [, dprcat=dprcat] [, xoffse=xoffse] [, yoffse=yoffse] ; [, angles=angles] ; ; REQUIRED INPUTS: ; fileli list of MIDI files ; keywrd keyword(s) in the FITS header of the MIDI files ; ; OPTIONAL INPUTS: ; ytitle string containing a user specified title for the y axis ; ; KEYWORDS: ; yrange range of values to be plotted on the y axis ; psfile name of the PS file for PS output ; titles title string for the entire plot ; labels string array containing user supplied lables ; dprcat string array containing the user supplied DPR category ; xoffse x offsets for the label positions ; xoffse y offsets for the label positions ; angles angles for plotting the label strings ; ; OUTPUTS: ; a plot of the values in KEYWRD with time of the observation ; ; DESCRIPTION AND EXAMPLE: ; This procedure makes a plot of the values stored in the FITS header ; keyword specified by the parameter KEYWRD for all the files given in ; the file list FILELI. If KEYWRD has more than one element, then the ; second element is interpreted as the error of the first element and ; error bars are also plotted. If the file list is of compressed MIDI ; type (i.e. several file names seprarated by blanks) then the first ; and last filename will be used to plot the specific keyword values. ; The title string for the entire plot (TITLES) and for the y axis ; (YTITLE), the position offsets (XOFFSE and YOFFSE) and rotation ; angles (ANGLES) for the label of each data point as well as the ; label strings (LABLES) can be specified explicitely by the caller. ; Common keywords to be plotted are for example: ; midi_pltprp, fileli, 'ISS PBLA12 END', 'PA ['+string("232B)+']' ; midi_pltprp, fileli, 'ISS PBL12 END', 'BL [m]' ; midi_pltprp, fileli, 'AIRM START', 'airmass', yrange=[1.6, 1.0] ; midi_pltprp, fileli, 'AMBI FWHM END', 'seeing [arcsec]' ; midi_pltprp, fileli, 'AMBI TAU0 END', '!4!X!D0!N [sec]' ; midi_pltprp, fileli, ['AO2 STREHL_MEAN', 'AO2 STREHL_RMS'], 'strehl' ; ; CALLED BY: ; midi_pltpds ; ; CALLING: ; setplot ; ; MODIFICATION HISTORY: ; 2008-10-14 Written by Konrad R. Tristram ; 2009-10-09 Konrad R. Tristram: Added possibility to plot errors. ; PRO MIDI_PLTPRP, FILELI, KEYWRD, YTITLE, PSFILE=PSFILE, YRANGE=YRANGE, $ TITLES=TITLES, LABELS=LABELS, DPRCAT=DPRCAT, XOFFSE=XOFFSE, $ YOFFSE=YOFFSE, ANGLES=ANGLES ; CHECK IF AT LEAST TWO PARAMETERS HAVE BEEN PASSED TO THE FUNCTION ;------------------------------------------------------------------------------- if n_params() lt 2 then begin message, 'Syntax - MIDI_PLTPRP, fileli, keywrd [, ytitle, ' + $ 'psfile=psfile,', /continue message, ' titles=titles, labels=labels,', /contin message, ' dprcat=dprcat, xoffse=xoffse,', /contin message, ' yoffse=yoffse, angles=angles]', /contin return endif ; CHECK AND INITIALISE SEVERAL PARAMETERS ;------------------------------------------------------------------------------- if (n_elements(ytitle) ne 1) then ytitle = strlowcase(keywrd) if (n_elements(xoffse) ne n_elements(fileli)) then $ xoffse = replicate(0.0075, n_elements(fileli)) if (n_elements(yoffse) ne n_elements(fileli)) then $ yoffse = replicate(-0.005, n_elements(fileli)) if (n_elements(angles) ne n_elements(fileli)) then $ angles = replicate(0.0000, n_elements(fileli)) keytmp = ['OBS TARG NAME', 'DPR CATG', 'DET NRTS MODE', 'MJD-OBS', keywrd] ; INITIALISE THE LABEL_DATE FUNCTION ;------------------------------------------------------------------------------- tmpvar = label_date(date_format='%H:%I') ; SELECT THE FIRST AND LAST FILENAME FOR EACH ELEMENT IN THE FILE LIST ;------------------------------------------------------------------------------- filbeg = fileli & filend = fileli for i=0, n_elements(fileli)-1 do begin $ tmpvar = strsplit(fileli[i], /extract, count=count) &$ filbeg[i] = tmpvar[0] &$ filend[i] = tmpvar[count-1] &$ endfor ; GET THE HEADER KEYWORDS FOR THE TWO FILE LISTS ;------------------------------------------------------------------------------- filestru = obj_new('fitsfilelist', filbeg) keybeg = filestru->getheaderkeys(keytmp) obj_destroy, filestru filestru = obj_new('fitsfilelist', filend) keyend = filestru->getheaderkeys(keytmp) obj_destroy, filestru ; SET THE LABELS FOR THE DATA POINTS IF NOT SUPPLIED ;------------------------------------------------------------------------------- if (n_elements(labels) ne n_elements(fileli)) then labels = keyend.obstargname if (n_elements(dprcat) ne n_elements(fileli)) then dprcat = keyend.dprcatg ; FIDDLE WITH THE DATA, E.G. WITH THE POSITION ANGLES ;----------------------------------------------------------------------- if keytmp[4] eq 'ISS PBLA12 END' then begin selidx = where(keybeg.(5) gt 210) if selidx[0] ge 0 then begin keybeg[selidx].(5) -=360 keyend[selidx].(5) -=360 endif endif ; CALCULATE THE X AXIS RANGE TO BE PLOTTED ;------------------------------------------------------------------------------- mindate = min(keybeg.mjdobs) + 2400000.5D maxdate = max(keyend.mjdobs) + 2400000.5D tmpdate = mindate mindate = 1.05D * tmpdate - 0.05D * maxdate maxdate = 1.05D * maxdate - 0.05D * tmpdate mindate = tmpdate ; CALCULATE THE Y AXIS RANGE TO BE PLOTTED ;------------------------------------------------------------------------------- if n_elements(yrange) ne 2 then begin minvalu = min([keybeg.(5), keyend.(5)]) maxvalu = max([keybeg.(5), keyend.(5)]) tmpvalu = minvalu minvalu = 1.05 * tmpvalu - 0.05 * maxvalu maxvalu = 1.05 * maxvalu - 0.05 * tmpvalu yrange = [minvalu, maxvalu] endif ; SET THE PLOTTING ENVIRONMENT, MAKE USER SYMBOL AND LOAD THE COLOUR TABLE ;------------------------------------------------------------------------------- setplot, pltset, pltpos, /normal, psfile=psfile usersym, cos(findgen(17) * (!pi*2/16.)), sin(findgen(17) * (!pi*2/16.)), /FILL if keyword_set(psfile) then symsiz=0.5 else symsiz=1.0 loadct, 39, /silent ; PLOT THE AXIS, THE LABLES AND TITLES ;------------------------------------------------------------------------------- plot, [mindate, maxdate], yrange, yrange=yrange, $ xticku = ['Hours'], xminor = 6, xtickformat='label_date', $ xtitle = 'time of observation [hrs:min]', ytitle=ytitle, $ title = titles, position=pltpos, /nodata, /device, /noerase xrange = !x.crange[1] - !x.crange[0] yrange = !y.crange[1] - !y.crange[0] ; LOOP OVER ALL ELEMENTS TO BE PLOTTED ;------------------------------------------------------------------------------- for i=0,n_elements(keybeg)-1 do begin ; SET THE COLOR DEPENDING ON THE DATA CATEGORY ;----------------------------------------------------------------------- if dprcat[i] eq 'SCIENCE ' then colour=60 else colour=150 ; SET THE SYMBOL DEPENDING ON THE DATA TYPE ;----------------------------------------------------------------------- case keybeg[i].detnrtsmode of 'ACQ_UT_COARSE_CHOP': symbol=-7 'OBS_FRINGE_SEARCH_FOURIER': symbol=-2 'OBS_FRINGE_SEARCH_DISPERSED': symbol=-2 'OBS_FRINGE_TRACK_FOURIER': symbol=-8 'OBS_FRINGE_TRACK_DISPERSED': symbol=-8 'OBS_FRINGE_TRACK_DISPERSED_OFF': symbol=-8 'DEFAULT_CHOP': symbol=-6 else: symbol=-5 endcase ; PLOT THE DATA ;----------------------------------------------------------------------- oplot, [keybeg[i].mjdobs, keyend[i].mjdobs] + 2400000.5D, $ [keybeg[i].(5), keyend[i].(5)], color=colour, psym = symbol, $ symsize = symsiz, thick=2 &$ ; IF SECOND KEYWORD PRESENT, INTERPRET IT AS ERRORS AND PLOT THESE ;----------------------------------------------------------------------- if n_tags(keybeg) eq 7 then begin errplot, [keybeg[i].mjdobs, keyend[i].mjdobs] + 2400000.5D, $ [keybeg[i].(5)-keybeg[i].(6), $ keyend[i].(5)-keyend[i].(6)], $ [keybeg[i].(5)+keybeg[i].(6), $ keyend[i].(5)+keyend[i].(6)], color=colour endif ; ADD A LABEL TO THE DATA POINT ;----------------------------------------------------------------------- xyouts, keyend[i].mjdobs + 2400000.5D + xoffse[i] * xrange, $ keyend[i].(5) + yoffse[i] * yrange, labels[i], $ charsize=0.75*symsiz, orientation=angles[i], $ color=colour endfor ; CALCULATE THE POSITIONS FOR A LEGEND ;------------------------------------------------------------------------------- xoffse = 0.05*(!x.crange[1]-!x.crange[0]) xposit = !x.crange[0] + 1.0 * xoffse yoffse = 0.05*(!y.crange[1]-!y.crange[0]) yposit = !y.crange[0] + 19.0 * yoffse ; PLOT A LEGEND ;------------------------------------------------------------------------------- plots, xposit+0.0*xoffse, yposit-0.0*yoffse, color=60, psym=8, $ symsize=symsiz, thick=2 xyouts, xposit+0.3*xoffse, yposit-0.2*yoffse, 'science' plots, xposit+0.0*xoffse, yposit-1.0*yoffse, color=150, psym=8, $ symsize=symsiz, thick=2 xyouts, xposit+0.3*xoffse, yposit-1.2*yoffse, 'calibrator' ; CLOSE THE PLOT ;------------------------------------------------------------------------------- setplot, pltset, /close END ;+ ; NAME: ; MIDI_PLTPDS ; ; PURPOSE: ; Plot common FITS header keywords for fringe tracks in a MIDI data set. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; midi_pltpds, datset, titles [, pspath=pspath] ; ; REQUIRED INPUTS: ; datset MIDI data set containing the fringe tracks ; ; OPTIONAL INPUTS: ; titles title string for the entire series of plots ; ; KEYWORDS: ; pspath prefix to the PS file names for PS output ; ; OUTPUTS: ; plots of several quantities with time of the observation ; ; DESCRIPTION AND EXAMPLE: ; This procedure is a wrapper for calling MIDI_PLTPRP for a MIDI data ; set. From the fringe track files given in the data set several plots ; are automatically created for the following quantities: projected ; baseline angle, projected baseline length, airmass, DIMM seeing and ; coherence time. If possible the quantities from the end of the file ; are used, e.g. the header keyword 'AIRM END' is used instead of ; 'AIRM START'. The procedure is called by for example: ; midi_pltpds, datset, 'My title' ; ; CALLED BY: ; none ; ; CALLING: ; midi_pltprp ; ; MODIFICATION HISTORY: ; 2008-10-14 Written by Konrad R. Tristram ; PRO MIDI_PLTPDS, DATSET, TITLES, PSPATH=PSPATH ; CHECK IF PS OUTPUT IS DESIRED AND SET VARIABLES ACCORDINGLY ;------------------------------------------------------------------------------- if keyword_set(pspath) then begin psuffx = ['_pangle', '_blength', '_airmass', '_seeing', '_tau0'] endif else begin pspath='' psuffx = replicate('', 5) endelse ; INITIALISE THE POSITIONING VARIABLES ;------------------------------------------------------------------------------- xoffse = replicate(-0.01, 6) yoffse = replicate(-0.02, 6) angles = replicate(-90.0, 6) ; PLOT THE PROJECTED BASELINE ANGLE ;------------------------------------------------------------------------------- midi_pltprp, datset.filnam[0], 'ISS PBLA12 END', 'projected baseline ' + $ 'angle PA ['+string("232B)+']', titles=titles, $ labels=datset.object, dprcat=datset.dprcat, xoffse=xoffse, $ yoffse=yoffse, angles=angles, psfile=pspath+psuffx[0] ; PLOT THE PROJECTED BASELINE LENGTH ;------------------------------------------------------------------------------- midi_pltprp, datset.filnam[0], 'ISS PBL12 END', 'projected baseline ' + $ 'length BL [m]', titles=titles, $ labels=datset.object, dprcat=datset.dprcat, xoffse=xoffse, $ yoffse=yoffse, angles=angles, psfile=pspath+psuffx[1] ; PLOT THE AIRMASS ;------------------------------------------------------------------------------- midi_pltprp, datset.filnam[0], 'ISS AIRM END', 'airmass', titles=titles, $ labels=datset.object, dprcat=datset.dprcat, xoffse=xoffse, $ yoffse=yoffse, angles=angles, psfile=pspath+psuffx[2] ; PLOT THE DIMM SEEING ;------------------------------------------------------------------------------- midi_pltprp, datset.filnam, 'ISS AMBI FWHM END', 'DIMM seeing ' + $ '[arcsec]', yrange=[0,2.0], titles=titles, $ labels=datset.object, dprcat=datset.dprcat, xoffse=xoffse, $ yoffse=yoffse, angles=angles, psfile=pspath+psuffx[3] ; PLOT THE COHERENCE TIME ;------------------------------------------------------------------------------- midi_pltprp, datset.filnam, 'ISS AMBI TAU0 END', 'coherence time ' + $ '!4s!X!D0!N [sec]', titles=titles, $ labels=datset.object, dprcat=datset.dprcat, xoffse=xoffse, $ yoffse=yoffse, angles=angles, psfile=pspath+psuffx[4] END ;+ ; NAME: ; MIDI_A2JAVE ; ; PURPOSE: ; Calculate the average and scatter of conversion factors. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; a2jave = midi_a2jave(adu2jy [, mjdobs] [, nummax] [, sigclp=sigclp] ; [, flgdat=flgdat] [,/silent]) ; ; REQUIRED INPUTS: ; adu2jy structure containing the ADU to Jy conversion factors ; ; OPTIONAL INPUTS: ; mjdobs a Julian date around which to select calibrators ; nummax maximum number of calibrators to select for averaging ; ; KEYWORDS: ; sigclp number of sigmas for clipping of calibrators ; silent suppress output of messages ; caldatabase IDL object containing the calibrator database ; silent suppress the output of infromative messages ; ; OUTPUTS: ; adu2jy the average and scatter of the conversion factors ; flgdat flagging information for the calibrators ; ; DESCRIPTION AND EXAMPLE: ; This function uses a set of conversion factors for calibrators (as ; for example created by the function MIDI_ADU2JY) to calculate the ; mean conversion factor and its scatter. If no Julian date is pro- ; vided through MJDOBS, then the calculation is carried out for all ; calibrators of the input structure. Otherwise, only the NUMMAX ; closest calibrators in time to the Julian date are selected for the ; calculation. Additionally all calibrators from different nights are ; discarded unless this would mean that less than MINCAL=2 calibrators ; remain. The selection is carried out independently for the photome- ; tries, the correlated flux, the visibility and the phase. Then a ; simple average and the standard deviation of the conversion factors ; is calculated. Iteratively, outliers of more than SIGCLP sigma are ; rejected and the new average and scatter calculated. The resulting ; values are returned by the function. In addition, the flagging data ; can be retrieved as a structure similar to ADU2JY by specifying the ; keyword FLGDAT. For every element of the input, and for each of the ; 5 quantities (masked photometry, unmasked photometry, correlated ; flux, the visibility and the phase) the flag is 0 if the data was ; used for the calculation, 1 if the data was rejected because of too ; large distance in time or incomplete data and 2 if the data was re- ; jected because it is an outlier. ; The function is easiest called by: ; a2jave = midi_a2jave(adu2jy) ; ; CALLED BY: ; midi_adu2jy ; midi_pltres ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2012-06-12 Written by Konrad R. Tristram ; 2012-08-29 Konrad R. Tristram: Use also information from QUALIT tag. ; FUNCTION MIDI_A2JAVE, ADU2JY, MJDOBS, NUMMAX, SIGCLP=SIGCLP, FLGDAT=FLGDAT, $ SILENT=SILENT ; SET DEFAULT VALUES OF A FEW VARIABLES ;------------------------------------------------------------------------------- if n_elements(nummax) lt 1 then nummax = 7 else nummax = nummax[0] if ~ keyword_set(silent) then silent = 0 mincal = 2 ; minimum allowed number of calibrators. ; GET THE INDICES OF THE TAGS THAT NEED TO BE PROCESSED ;------------------------------------------------------------------------------- tagnam = ['MSKFLX', 'MSKFLXERR', 'TOTFLX', 'TOTFLXERR', 'CORFLX', 'CORFLXERR', $ 'VISAMP', 'VISAMPERR', 'VISPHI', 'VISPHIERR'] tagidx = intarr(n_elements(tagnam)) for i=0,n_elements(tagnam)-1 do tagidx[i]=where(tag_names(adu2jy) eq tagnam[i]) ; CREATE THE STRUCTURE HOLDING THE FLAGS ;------------------------------------------------------------------------------- flgdat = {identi:'', object:''} for i=0,n_elements(tagnam)/2-1 do flgdat=create_struct(flgdat, tagnam[2*i], 1B) flgdat = replicate(flgdat, n_elements(adu2jy)) for i=0,1 do flgdat.(i) = adu2jy.(i) ; CREATE THE STRUCTURE HOLDING THE AVERAGE AND ERROR OF THE CONVERSION FACTORS ;------------------------------------------------------------------------------- a2jave = adu2jy[0] for i=0,n_elements(tag_names(a2jave))-1 do a2jave.(i) = '' a2jave.wavele = adu2jy[0].wavele a2jave.object = 'Variation of conversion factors' ; DECOMPOSE THE QUALITY FLAGS INTO AN ARRAY ;------------------------------------------------------------------------------- qualit = reverse(byte((string(adu2jy.qualit, format='(B08)')))-48B) qualit = qualit[[2,2,1,3,1], *] ; LOOP OVER THE DIFFERENT TAGS ;------------------------------------------------------------------------------- for i=0,4 do begin ; PRINT A DIAGNOSTIC MESSAGE ;------------------------------------------------------------------------------- if ~ silent then message, '### Processing the '+tagnam[2*i]+' data.', /continue ; CHECK WHICH CALIBRATORS ACTUALLY CONTAIN USEFUL DATA AND ARE UNFLAGGED ;------------------------------------------------------------------------------- wavsel = where(adu2jy[0].wavele gt 8.2 and adu2jy[0].wavele lt 13.0) select = where((total(adu2jy.(tagidx[2*i])[wavsel], 1) ne 0) and $ (qualit[i, *] eq 0)) ; CHECK IF THERE ARE ENOUGH ELEMENTS TO CALCULATE A VARIATION ;------------------------------------------------------------------------------- if n_elements(select) ge mincal then begin ; CHECK IF A REFERENCE DATE WAS PROVIDED ;----------------------------------------------------------------------- if n_elements(mjdobs) gt 0 then begin ; FIND THE NUMMAX CLOSEST CALIBRATORS IN TIME TO MJDOBS ;--------------------------------------------------------------- tmpvar = abs(adu2jy[select].mjdobs - mjdobs) sorted = sort(tmpvar) sorted = sorted[0:(n_elements(sorted)-1<(nummax-1)>(mincal-1))] ; CHECK IF THERE ARE TIME DIFFERENCES OF MORE THAN 1 NIGHT ;--------------------------------------------------------------- badcal = (tmpvar[sorted] gt 0.6) badidx = where(badcal) ; PRINT INFORMATION OR WARNING IF NOT SUPPRESSED ;--------------------------------------------------------------- if ~ silent then begin if badidx[0] ge 0 then for j=0,n_elements(badidx)-1 do $ message, string(select[sorted[badidx[j]]], $ adu2jy[select[sorted[badidx[j]]]].object, $ adu2jy[select[sorted[badidx[j]]]].mjdobs, $ format='("Calibrator", I3, ": ", A12, " at ' + $ 'MJD =", F12.5, " is from a different ' + $ 'night! Ignoring it.")'), /continue if (total(badcal[0:mincal-1]) gt 0) then $ message, string(mincal, format='("Warning: Found less ' + $ 'than", I2, " calibrators in the same night! ' + $ 'Retaning the two closest in time.")'), /continue endif ; SELECT ONLY THE GOOD CALIBRATORS ;--------------------------------------------------------------- badcal[0:mincal-1] = 0 select = select[sorted(where(~ badcal))] endif flgdat[select].(i+2) = 0. ; SET DEFAULT VALUE FOR OUTLIER REJECTION ;----------------------------------------------------------------------- if n_elements(sigclp) lt 1 then begin if n_elements(select) ge 5 then numsig = 2.0 else numsig = 1.5 endif else numsig = sigclp[0] > 0.1 ; LOOP OVER OUTLIER REJECTION ;----------------------------------------------------------------------- badidx = 0 while badidx[0] ge 0 do begin ; CALCULATE THE MEAN TRANSFER FUNCTION BY SIMPLY AVERAGING ;--------------------------------------------------------------- averag = total(adu2jy[select].(tagidx[2*i]), 2)/$ n_elements(select) ; CALCULATE THE SCATTER OF THE TRANSFER FUNCTIONS ;--------------------------------------------------------------- deviat = fltarr(n_elements(a2jave.wavele)) for j=0,n_elements(a2jave.wavele)-1 do begin deviat[j] = stddev(adu2jy[select].(tagidx[2*i])[j]) endfor ; FIND OUTLIERS BY COMPARING THE DEVIATIONS TO THE SCATTER ;--------------------------------------------------------------- badcal = averag[wavsel] # replicate(1, n_elements(select)) badcal = (adu2jy[select].(tagidx[2*i])[wavsel] - badcal)^2 badcal = total(sqrt(badcal), 1) gt numsig*total(deviat[wavsel]) badidx = where(badcal) ; badidx = where(badcal and ~(qualit[i, select])) ; this is the old code, but seems erroneous ; PRINT A MESSAGE AND UPDATE FLAGGING IF OUTLIERS WERE FOUND ;--------------------------------------------------------------- if badidx[0] ge 0 then for j=0,n_elements(badidx)-1 do begin if (~ silent) then message, string(select[badidx[j]], $ adu2jy[select[badidx[j]]].object, $ adu2jy[select[badidx[j]]].mjdobs, numsig, $ format='("Calibrator", I3, ": ", A12, " at ' + $ 'MJD =", F12.5, " deviates by more than", ' + $ 'F4.1, " sigma! Ignoring it.")'), /continue flgdat[select[badidx[j]]].(i+2) = 2. endfor ; SELECT ONLY THE REMAINING GOOD CALIBRATORS ;--------------------------------------------------------------- select = select[where(~ badcal)] endwhile ; COPY THE AVERAGE AND THE SCATTER INTO THE RESULT STRUCTURE ;----------------------------------------------------------------------- a2jave.(tagidx[2*i ]) = averag a2jave.(tagidx[2*i+1]) = deviat ; IF THERE ARE NOT ENOUGH ELEMENTS SET MEAN AND ERROR TO ZERO ;------------------------------------------------------------------------------- endif else begin a2jave.(tagidx[2*i ]) = 0. a2jave.(tagidx[2*i+1]) = 0. endelse ; END OF LOOP OVER THE DIFFERENT TAGS ;------------------------------------------------------------------------------- endfor ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, a2jave END ;+ ; NAME: ; MIDI_ADU2JY ; ; PURPOSE: ; Calculate the ADU to JY conversion factor for a set of calibrators. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; adu2jy = midi_adu2jy(datset [, titles] [, pspath=pspath] [,/noplot] ; [, caldatabase=caldatabase] [,/silent]) ; ; REQUIRED INPUTS: ; datset MIDI data set containing the EWS tags for fringe tracks ; ; OPTIONAL INPUTS: ; titles title string for the entire series of plots ; ; KEYWORDS: ; pspath prefix to the PS file names for PS output ; noplot suppress the output of plots ; plterr also plot error bars for the data plotted ; caldatabase IDL object containing the calibrator database ; sigclp number of sigmas for clipping of calibrators ; silent suppress the output of infromative messages ; ; OUTPUTS: ; three plots of the ADU to JY conversion factor ; adu2jy struct containing the ADU to Jy conversion factor ; ; DESCRIPTION AND EXAMPLE: ; This function loads the raw data corresponding to the MIDI data set ; given in DATSET. For each object in the raw data / data set, a ; matching calibrator in a calibrator database is determined by ; comparing the coordinates of the FITS file to those of the entries ; in the calibrator database. By default the calibrator database by ; Roy van Boekel (built into MIA+EWS) is used. If a match is found, ; then the spectrum in the database is used as the template spectrum. ; Otherwise a Rayleigh-Jeans law, scaled using the calibrator flux ; specified in DATSET, is used as the template spectrum. If neither a ; match in the database is found nor a calibrator flux is given in ; DATSET, then the object is skipped. The raw total flux, the raw ; masked flux and the raw correlated flux are then devided by the ; template to obtain the respective conversion factor. The average ; and the scatter of the conversion factors are calculated by calling ; the function MIDI_A2JAVE. Finally three plots are created with the ; conversion factors for all of the objects as well as their average ; (dashed) and the scatter (dotted). The function returns a struct ; containing the conversion factors for the individual calibrators as ; well as the mean and standard deviation of the conversion factors. ; In order to run the programme in batch mode it is possible to ; suppress the creation of plots using the keyword NOPLOT and to ; provide the calibrator database through the keyword CALDATABASE. ; The function may be called by: ; adu2jy = midi_adu2jy(datset) ; ; CALLED BY: ; none ; ; CALLING: ; midi_getdat ; midi_a2jave ; setplot ; ; MODIFICATION HISTORY: ; 2008-10-14 Written by Konrad R. Tristram ; 2008-11-06 Konrad R. Tristram: Changed text output of the programme. ; 2009-05-11 Konrad R. Tristram: Added plot of standard deviation and ; removal of '0' before calibrator number. ; 2009-05-12 Konrad R. Tristram: Modified to cope with new CALCAT struct. ; 2010-07-20 Konrad R. Tristram: Modified to use calib database in EWS. ; 2010-12-13 Konrad R. Tristram: Add MJDOBS in the resulting structure. ; 2012-04-20 Konrad R. Tristram: Also calculate and return variation. ; 2012-04-23 Konrad R. Tristram: Added CALDATABASE, NOPLOT, SILENT. ; 2012-06-12 Konrad R. Tristram: Rewrite, changed procedure to a function. ; 2012-08-18 Konrad R. Tristram: Modified labelling of the calibrators. ; 2012-08-29 Konrad R. Tristram: Use also information from QUALIT tag. ; 2012-10-11 Konrad R. Tristram: Added keyword SIGCLP for MIDI_A2JAVE. ; FUNCTION MIDI_ADU2JY, DATSET, TITLES, PSPATH=PSPATH, NOPLOT=NOPLOT, $ PLTERR=PLTERR, CALDATABASE=CALDATABASE, $ SIGCLP=SIGCLP, SILENT=SILENT ; INITIALISE A FEW VARIABLES ;------------------------------------------------------------------------------- tagnam = ['MSKFLX', 'MSKFLXERR', 'TOTFLX', 'TOTFLXERR', 'CORFLX', 'CORFLXERR'] tprefx = ['Masked flux', 'Total flux', 'Correlated flux']+' conversion factor' ytitle = ['ADU / Jy (masked flux)', 'ADU / Jy (total flux)', $ 'ADU / Jy (correlated flux)'] if ~ keyword_set(silent) then silent = 0 if ~ keyword_set(plterr) then plterr = 0 ; CHECK IF AN ADDITIONAL TITLE STRING IS PROVIDED AND CREATE THE FINAL TITLES ;------------------------------------------------------------------------------- if keyword_set(titles) then titstr = tprefx + ': ' + titles $ else titstr = tprefx ; CHECK IF PS OUTPUT IS DESIRED AND SET VARIABLES ACCORDINGLY ;------------------------------------------------------------------------------- if keyword_set(pspath) then psuffx = ['_adu2msk', '_adu2tot', '_adu2cor'] $ else begin psuffx = ['', '', ''] & pspath = '' & endelse ; READ IN THE RAW DATA ;------------------------------------------------------------------------------- rawdat = midi_getdat(datset, /raw) ; CREATE THE STRUCTURE HOLDING THE RESULTING CONVERSION FACTORS ;------------------------------------------------------------------------------- adu2jy = rawdat tagidx = intarr(n_elements(tagnam)) for i=0,5 do begin tagidx[i] = where(tag_names(adu2jy) eq tagnam[i]) if tagidx[i] ge 0 then adu2jy.(tagidx[i]) = 0.0 endfor ; LOAD THE CALIBRATOR DATA BASE OF ROY VAN BOEKEL IF NOT PROVIDED ;------------------------------------------------------------------------------- if size(caldatabase, /type) ne 11 then begin boecal = vboekelbase() endif else begin boecal = caldatabase endelse ; LOOP OVER ALL ELEMENTS IN THE RAW DATA ;------------------------------------------------------------------------------- for i=0, n_elements(rawdat)-1 do begin ; DECOMPOSE THE QUALITY FLAGS INTO AN ARRAY ;----------------------------------------------------------------------- qualit = reverse(byte((string(adu2jy[i].qualit, format='(B08)')))-48B) ; FIND CALIBRATOR WITHIN 2 ARCMIN FROM POSITION IN FITS FILE ;----------------------------------------------------------------------- distan = boecal->sourcesNearPosition(adu2jy[i].rascen, $ adu2jy[i].declin, srcdat, nmax=1, rmax=1./30., ierr=errvar) ; IF A CALIBRATOR WAS FOUND PRINT A MESSAGE AND INTERPOLATE THE FLUX ;----------------------------------------------------------------------- if not errvar then begin tmpvar = boecal->sourceByName(srcdat[0].name, srcdat, diamet, $ photom, spectr, ierr=errvar) if ~ silent then print, i, rawdat[i].object, $ strtrim(srcdat.name, 2), format='(I2, ": ", A10, ' + $ '" found in the calibrator data base as ", A, ". ' + $ 'Using the modelled spectrum.")' spectr = *(spectr[0]) templa = interpol(spectr.flux, 1.E6 * spectr.wavelength, $ rawdat[i].wavele) qualit[4] = 1B tplerr = 0.05 ; estimated error of the calibrator data base ; IF NO CALIBRATOR WAS FOUND BUT HAVE A FLUX, MAKE R.-J. APPROXIMTION ;----------------------------------------------------------------------- endif else begin if datset[i].calflx gt 0.0 then begin if ~ silent then print, i, rawdat[i].object, datset[i].calflx, $ format='(I2, ": ", A10, " not found in the ' + $ 'calibrator data base. Using Rayleigh-Jeans ' + $ 'approximation with F =", F6.2, " Jy.")' templa = rawdat[i].wavele^(-2) * 100. * datset[i].calflx qualit[4:5] = [0B, 1B] tplerr = 0.10 ; estimated error of the flux provided ; OTHERWISE WE HAVE NO TEMPLATE SO WE SET THE QUALITY FLAGS TO 0 ;----------------------------------------------------------------------- endif else begin if ~ silent then print, i, rawdat[i].object, $ format='(I2, ": ", A10, " not found in the ' + $ 'calibrator data base and no flux given. ' + $ 'Skipping the source.")' qualit[4:5] = [0B, 0B] endelse & endelse ; IF THERE IS A TEMPLATE CALCULATE THE CONVERSION FACTORS ;----------------------------------------------------------------------- if total(qualit[4:5]) gt 0 then begin adu2jy[i].totflx = rawdat[i].totflx / templa adu2jy[i].totflxerr = sqrt((rawdat[i].totflxerr/templa)^2 + $ (tplerr * adu2jy[i].totflx)^2) adu2jy[i].mskflx = rawdat[i].mskflx / templa adu2jy[i].mskflxerr = sqrt((rawdat[i].mskflxerr/templa)^2 + $ (tplerr * adu2jy[i].mskflx)^2) ; There should be a correction for partially resolved objects!! adu2jy[i].corflx = rawdat[i].corflx / templa adu2jy[i].corflxerr = sqrt((rawdat[i].corflxerr/templa)^2 + $ (tplerr * adu2jy[i].corflx)^2) endif ; COPY BACK THE QUALITY FLAGS ;----------------------------------------------------------------------- adu2jy[i].qualit = total(qualit * 2^(indgen(8))) endfor ; DESTROY THE CALIBRATOR DATABASE IF NOT COMING FROM OUTSIDE THE ROUTINE ;------------------------------------------------------------------------------- if size(caldatabase, /type) ne 11 then obj_destroy, boecal ; ONLY SELECT THE OBJECTS WHERE A CONVERSION FACTOR WAS DETERMINED ;------------------------------------------------------------------------------- select = byte(string(adu2jy.qualit, format='(B08)'))-48B select = where(total(select[2:3,*], 1) gt 0) if select[0] lt 0 then begin message, 'No calibrators found in the data. Exiting.', /continue return, -1 endif else adu2jy = adu2jy[select] numele = n_elements(adu2jy) ; DETERMINE THE AVERAGE AND SCATTER OF THE TRANSFER FUNCTION ;------------------------------------------------------------------------------- a2jave = midi_a2jave(adu2jy, flgdat=flgdat, sigclp=sigclp, /silent) ; ADD THE FLAGGING INFORMATION TO THE QUALITY TAG ;------------------------------------------------------------------------------- select = [2, 2, 1, 3] qualit = reverse(byte((string(adu2jy.qualit, format='(B08)')))-48B) for i=0,3 do begin badstr = where(flgdat.(i+2) ge 1) if badstr[0] ge 0 then qualit[select[i], badstr] = 1B endfor adu2jy.qualit = total(qualit*(2^(indgen(8))#replicate(1,n_elements(qualit))), 1) ; LOOP OVER THE THREE TYPES OF PLOTS ;------------------------------------------------------------------------------- if ~ keyword_set(noplot) then for i=0,2 do begin ; OPEN THE PLOT DEVICE AND CREATE THE PLOT ;----------------------------------------------------------------------- setplot, pltset, /normal, psfile=pspath + psuffx[i] yrange = [0, max(adu2jy.(tagidx[2*i]) + $ (adu2jy.(tagidx[2*i+1]) < (0.5*adu2jy.(tagidx[2*i]))))] plot, [5, 14], yrange, xstyle=5, ystyle=5, /nodata, /noerase loadct, 33, /silent ; DETERMINE THE INDEX OF THE CURRENT TAG IN THE FLAGGING STRUCT ;----------------------------------------------------------------------- flgidx = where(tag_names(flgdat) eq tagnam[2*i]) ; LOOP OVER THE INDIVIDUAL OBJECTS TO PLOT AND LABEL THEM ;----------------------------------------------------------------------- for j=0,numele-1 do begin ; PLOT THE DATA AND ITS ERRORS ;--------------------------------------------------------------- colour = 165 + (1.0*j/(numele-1))^1.2 * (254-164) oplot, adu2jy[j].wavele, adu2jy[j].(tagidx[2*i]), color=colour if plterr then begin errplot, adu2jy[j].wavele, color=colour, width=0.0, $ adu2jy[j].(tagidx[2*i])-adu2jy[j].(tagidx[2*i+1]), $ adu2jy[j].(tagidx[2*i])+adu2jy[j].(tagidx[2*i+1]) endif ; ADD A STAR TO R.-J. APPROXIMATIONS ;--------------------------------------------------------------- if qualit[5,j] then badstr = ['', '*'] $ else badstr = ['', ''] ; ADD BRACKETS TO REJECTED CALIBRATORS ;--------------------------------------------------------------- if flgdat[j].(flgidx) ge 1 then badstr += ['(', ')'] $ else badstr += [' ', ''] ; PRINT OUT THE CALIBRATOR NAME ;--------------------------------------------------------------- xyouts, 5.4, (0.93-j/20.)*!y.crange[1], badstr[0]+ $ strtrim(adu2jy[j].object,2)+badstr[1], color=colour endfor ; PLOT THE AVERAGE AND THE ERROR OF THE RESPECTIVE CONVERSION FACTOR ;----------------------------------------------------------------------- loadct, 39, /silent oplot, a2jave.wavele, a2jave.(tagidx[2*i ]), thick=2, linestyle=5 oplot, a2jave.wavele, a2jave.(tagidx[2*i+1]), thick=2, linestyle=1 ; OVERPLOT THE AXIS AND LABELS ;----------------------------------------------------------------------- plot, [5, 14], yrange, xstyle = 1, /nodata, ytitle=ytitle[i], $ xtitle='wavelength !4k!X [!4l!Xm]', title=titstr[i], /noerase ; CLOSE THE PLOT ;----------------------------------------------------------------------- setplot, pltset, /close endfor ; PREPARE THE VALUE TO BE RETURNED ;------------------------------------------------------------------------------- return, [adu2jy, a2jave] END ;+ ; NAME: ; MIDI_CALERR ; ; PURPOSE: ; Estimate and add the errors from the variation of the transfer function. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; newdat = midi_calerr(caldat, datset[, mincal=mincal][, maxcal=maxcal]) ; ; REQUIRED INPUTS: ; caldat MIDI result structure containing the calibrated science data ; datset MIDI data set structure containing the calibrator data ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; mincal minimum number of calibrators to use ; maxcal maximum number of calibrators to use ; ; OUTPUTS: ; calibs an array containing the identies of the calibrators used ; newcal struct containing the calibrated data with modified errors ; ; DESCRIPTION AND EXAMPLE: ; This function determines the variation of the transfer function at the ; time of a science observation and adds this variation to the error of ; the data in quadrature. In a first step, the transfer functions (ADU ; to Jy conversion factors) for all calibrators in DATSET are determined. ; Then, for each calibrated data set in CALDAT, a penalty function is ; calculated for these calibrators, depending on the time difference to ; the calibrated data set. The penalty function also takes into account ; incomplete calibrator data (i.e. calibrators without correlated flux ; or photometry). Then at least MINCAL and at most MAXCAL calibrators ; with the lowest penalty function are selected. From these, the vari- ; ation of the transfer function is calculated for the masked flux, the ; total flux, the correlated flux, the visibility and the differential ; phases. This variation is then added to the error of the respective ; quantity in quadrature. ; The function may be called by: ; newcal = midi_calerr(caldat, datset) ; ; CALLED BY: ; none ; ; CALLING: ; midi_adu2jy ; ; MODIFICATION HISTORY: ; 2012-06-12 Written by Konrad R. Tristram ; 2012-10-10 Konrad R. Tristram: Fixed bugs and finalised documentation. ; FUNCTION MIDI_CALERR, CALDAT, DATSET, CALIBS, MINCAL=MINCAL, MAXCAL=MAXCAL ; SET THE DEFAULT VALUE FOR THE MAXIMUM NUMBER OF CALIBRATORS TO USE ;------------------------------------------------------------------------------- if n_elements(mincal) ne 1 then mincal = 3 if n_elements(maxcal) ne 1 then maxcal = 5 ; SELECT THE WAVELENGTH RANGE CONTAINING USEFUL DATA ;------------------------------------------------------------------------------- wavsel = where(caldat[0].wavele gt 8.0 and caldat[0].wavele lt 13.0) ; COMPILE A LIST OF ALL (USEFUL) CALIBRATORS IN DATSET ONLY ;------------------------------------------------------------------------------- adu2jy = midi_adu2jy(datset, /noplot, /silent) ; IF THERE ARE LESS THAN 3 CALIBRATORS CANNOT DO ANYTHING ;------------------------------------------------------------------------------- tmpvar = n_elements(adu2jy) if tmpvar lt 4 then begin message, 'Not enough useful calibrators found to calculate the ' + $ 'variation of the transfer function. Exiting!', /continue return, -1 endif else adu2jy = adu2jy[0:tmpvar-2] ; COPY THE INPUT DATA INTO THE RESULTING DATA ;------------------------------------------------------------------------------- newdat = caldat ; MAKE THE ARRAY CONTAINING THE IDENTITIES OF THE CALIBRATORS USED ;------------------------------------------------------------------------------- calibs = strarr(n_elements(caldat), maxcal) ; LOOP OVER ALL ELEMENTS IN THE CALIBRATED DATA ;------------------------------------------------------------------------------- for i=0, n_elements(caldat)-1 do begin ; CALCULATE PENALTY FOR TIME DIFFERENCE ;----------------------------------------------------------------------- penalt = abs(caldat[i].mjdobs - adu2jy.mjdobs) ; ADD PENALTY FOR MISMATCH IN COMPLETENESS OF CORRELATED FLUX ;----------------------------------------------------------------------- if total(caldat[i].corflx[wavsel], 1) ne 0 then begin select = total(adu2jy.corflx[wavsel], 1) eq 0 penalt += 1E6 * select endif ; ADD PENALTY FOR MISMATCH IN COMPLETENESS OF PHOTOMETRIES ;----------------------------------------------------------------------- if total(caldat[i].totflx[wavsel], 1) ne 0 then begin select = total(adu2jy.totflx[wavsel], 1) eq 0 penalt += 1E6 * select endif ; ADD PENALTY FOR DIFFERENT BASELINES ;----------------------------------------------------------------------- ; In principle more penalties should be added here, e.g. for different ; instrumental settings, baselines or airmass. ; DETERMINE THE MAXCAL CALIBRATORS WITH THE LOWEST PENALTY ;----------------------------------------------------------------------- sorted = (sort(penalt)) maxidx = ((where(penalt[sorted] gt 0.5))[0] - 1) > (mincal-1) sorted = sorted[0:(maxidx < (maxcal-1))] ; SAVE THE IDENTITIES OF THE SELECTED CALIBRATORS ;----------------------------------------------------------------------- calibs[i, 0:n_elements(sorted)-1] = adu2jy[sorted].identi ; GET THE NUMBER OF WAVELENGTH BINS AND CREATE A TEMPORARY VARIABLE ;----------------------------------------------------------------------- wavnum = n_elements(caldat[i].wavele) deviat = fltarr(wavnum) ; CALCULATE THE ERROR OF THE MASKED FLUX AND ADD IT IN QUADRATURE ;----------------------------------------------------------------------- if total(caldat[i].mskflx[wavsel], 1) ne 0 then begin averag = total(adu2jy[sorted].mskflx, 2)/n_elements(sorted) for j=0,wavnum-1 do deviat[j] = stddev(adu2jy[sorted].mskflx[j]) deviat = abs(deviat/averag) * newdat[i].mskflx newdat[i].mskflxerr = sqrt(newdat[i].mskflxerr^2 + deviat^2) endif ; CALCULATE THE ERROR OF THE TOTAL FLUX AND ADD IT IN QUADRATURE ;----------------------------------------------------------------------- if total(caldat[i].totflx[wavsel], 1) ne 0 then begin averag = total(adu2jy[sorted].totflx, 2)/n_elements(sorted) for j=0,wavnum-1 do deviat[j] = stddev(adu2jy[sorted].totflx[j]) deviat = abs(deviat/averag) * newdat[i].totflx newdat[i].totflxerr = sqrt(newdat[i].totflxerr^2 + deviat^2) endif ; CALCULATE THE ERROR OF THE CORRELATED FLUX AND ADD IT IN QUADRATURE ;----------------------------------------------------------------------- if total(caldat[i].corflx[wavsel], 1) ne 0 then begin averag = total(adu2jy[sorted].corflx, 2)/n_elements(sorted) for j=0,wavnum-1 do deviat[j] = stddev(adu2jy[sorted].corflx[j]) deviat = abs(deviat/averag) * newdat[i].corflx newdat[i].corflxerr = sqrt(newdat[i].corflxerr^2 + deviat^2) endif ; CALCULATE THE ERROR OF THE VISIBILITY AND ADD IT IN QUADRATURE ;----------------------------------------------------------------------- if total(caldat[i].visamp[wavsel], 1) ne 0 then begin averag = total(adu2jy[sorted].visamp, 2)/n_elements(sorted) for j=0,wavnum-1 do deviat[j] = stddev(adu2jy[sorted].visamp[j]) deviat = abs(deviat/averag) * newdat[i].visamp newdat[i].visamperr = sqrt(newdat[i].visamperr^2 + deviat^2) endif ; CALCULATE THE ERROR OF THE PHASES AND ADD IT IN QUADRATURE ;----------------------------------------------------------------------- if total(caldat[i].visphi[wavsel], 1) ne 0 then begin averag = total(adu2jy[sorted].visphi, 2)/n_elements(sorted) for j=0,wavnum-1 do deviat[j] = stddev(adu2jy[sorted].visphi[j]) ;deviat = abs(deviat/averag) * newdat[i].visphi newdat[i].visphierr = sqrt(newdat[i].visphierr^2 + deviat^2) endif endfor ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, newdat END ;+ ; NAME: ; MIDI_STRC2T ; ; PURPOSE: ; Print out an IDL structure as an ASCII table. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; midi_strc2t, datset [, tagnam] [, maxwid=maxwid] ; ; REQUIRED INPUTS: ; datset a MIDI data set or list of MIDI data sets ; ; OPTIONAL INPUTS: ; tagnam the names of the tags to be printed out as columns ; ; KEYWORDS: ; maxwid maximum with of a (text) column ; ; OUTPUTS: ; a table on screen ; ; DESCRIPTION AND EXAMPLE: ; This procedure prints out the contents of a struct array in form of a ; table of fixed size, with each column corresponding to one tag in the ; struct. The size of the column is determined from the data type and ; the space needed by the largest element in this column. For strings, ; the width of the column is limited by the keyword MAXWID; its default ; value is 40. The names of the tags and hence columns have to be ; specified by the string array TAGNAM. If TAGNAM is not set then all ; tags of the struct are printed as columns. For MIDI structures of ; type DATSET and RAWDAT/CALDAT a number of columns are automatically ; selected. The program is called by: ; midi_strc2t, datset ; ; CALLED BY: ; None. ; ; CALLING: ; strctxt ; ; MODIFICATION HISTORY: ; 2008-10-07 Written by Konrad R. Tristram ; 2008-11-06 Konrad R. Tristram: Added handling of BYTE type and arrays. ; 2009-06-26 Konrad R. Tristram: Rewritten using the programme strctxt. ; PRO MIDI_STRC2T, DATSET, TAGNAM, MAXWID=MAXWID ; GET THE TAG NAMES IN THE DATSET ;------------------------------------------------------------------------------- dattag = tag_names(datset) ; SELECT SPECIAL TAG NAMES FOR A FEW SPECIAL TYPES OF DATSET ;------------------------------------------------------------------------------- if n_elements(tagnam) eq 0 then begin if total('DPRCAT' eq dattag) eq 1 then begin tagnam=['IDENTI', 'FILNAM[0]', 'FILNAM[1]', 'FILNAM[2]', $ 'EWSTAG', 'DPRCAT', 'OBJECT', 'REDPAR.MSKFIL', $ 'REDPAR.TSMOOT', 'REDPAR.GSMOOT', 'REDPAR.DAVERA', $ 'REDPAR.SKYDIS', 'REDPAR.NOPHOT', 'REDPAR.CALIBR', $ 'CALFLX', 'CALDIA', 'COMMEN'] tagnam=['IDENTI', 'EWSTAG', 'DPRCAT', 'OBJECT', $ 'REDPAR.CALIBR', 'CALFLX', 'CALDIA', 'COMMEN'] endif else if total('UCOORD' eq dattag) eq 1 then begin tagnam=['IDENTI', 'OBJECT', 'BASELI', 'PBLINE', 'PANGLE', $ 'UCOORD', 'VCOORD', 'COMMEN'] endif else tagnam = dattag endif ; RUN STRCTXT TO PRINT OUT THE TABLE ;------------------------------------------------------------------------------- strctxt, datset, tagnam, maxwid=maxwid END ;+ ; NAME: ; MIDI_DATA2T ; ; PURPOSE: ; Print out MIDI data as an ASCII file ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; midi_data2t, caldat, filnam [, wavele] ; ; REQUIRED INPUTS: ; caldat a MIDI result structure containing the reduced data ; filnam the file name for the output file ; ; OPTIONAL INPUTS: ; wavele a two element array specifying the wavelength range ; ; KEYWORDS: ; none ; ; OUTPUTS: ; an ASCII file with the data ; ; DESCRIPTION AND EXAMPLE: ; This procedure prints out the contents of a MIDI result structure into ; an ASCII file. For every data set a small header including the object ; name and the uv coordinates is printed. Then the data is printed in ; columns. The program is called by: ; midi_data2t, caldat, filnam ; ; CALLED BY: ; None. ; ; CALLING: ; None. ; ; MODIFICATION HISTORY: ; 2009-06-23 Written by Konrad R. Tristram ; PRO MIDI_DATA2T, CALDAT, FILNAM, WAVELE ; SET THE DEFAULT VALUE FOR THE WAVELENGTH RANGE ;------------------------------------------------------------------------------- if n_elements(wavele) ne 2 then wavele = [7.5, 13.5] ; OPEN THE FILE FOR WRITING ;------------------------------------------------------------------------------- get_lun, lunnum openw, lunnum, filnam ; LOOP OVER THE UV POINTS IN THE DATA ;------------------------------------------------------------------------------- for i=0,n_elements(caldat)-1 do begin ; SELECT THE WAVELENGTHS IN THE WAVELENGTH RANGE ;----------------------------------------------------------------------- wavsel = where((caldat[i].wavele lt wavele[1]) and $ (caldat[i].wavele gt wavele[0])) ; PRINT THE HEADER FOR EVERY UV POINT ;----------------------------------------------------------------------- printf, lunnum, '# object : ' + caldat[i].object printf, lunnum, '# baseline : ' + caldat[i].baseli printf, lunnum, '# u coord [m]: ', caldat[i].ucoord printf, lunnum, '# v coord [m]: ', caldat[i].vcoord printf, lunnum, '# WAVELE TOTFLX TOTFLXERR CORFLX' + $ ' CORFLXERR VISAMP VISAMPERR VISPHI' + $ ' VISPHIERR' printf, lunnum, '# [micron] [Jy] [Jy] [Jy]' + $ ' [Jy] [degree]' + $ ' [degree]' ; PRINT THE DATA ;----------------------------------------------------------------------- for j=0,n_elements(wavsel)-1 do begin printf, lunnum, caldat[i].wavele[wavsel[j]], $ caldat[i].totflx[wavsel[j]], $ caldat[i].totflxerr[wavsel[j]], $ caldat[i].corflx[wavsel[j]], $ caldat[i].corflxerr[wavsel[j]], $ caldat[i].visamp[wavsel[j]], $ caldat[i].visamperr[wavsel[j]], $ caldat[i].visphi[wavsel[j]], $ caldat[i].visphierr[wavsel[j]], $ format='(9F12.6)' endfor printf, lunnum, '#' endfor ; CLOSE THE FILE ;------------------------------------------------------------------------------- close, lunnum END ;+ ; NAME: ; MIDI_DATMRG ; ; PURPOSE: ; Merge several data sets using Gaussian error propagation ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; result = midi_datmrg(caldat) ; ; REQUIRED INPUTS: ; caldat a MIDI result structure containing the reduced data ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; none ; ; OUTPUTS: ; result a result structure with the merged values ; ; DESCRIPTION AND EXAMPLE: ; This procedure calculates the average of several MIDI result datasets ; using Gaussian error propagation. As input, an array of structures ; with reduced MIDI data must be provided. This structure must contain ; certain keywords, such as CORFLX, MSKFLX, TOTFLX, VISAMP and VISPHI. ; These fields are averaged using weights determined from the statisti- ; cal errors, i.e.: ; result = total(value/error^2)/total(1/error^2) ; The errors of the result are calculated as follows: ; error = 1/sqrt(total(1/error^2)) ; The UV coordinates are averaged by determining the simple mean, the ; returned observation time is that of the first dataset. The function ; is called by: ; result = midi_datmrg(caldat) ; ; CALLED BY: ; None. ; ; CALLING: ; None. ; ; MODIFICATION HISTORY: ; 2011-01-26 Written by Konrad R. Tristram ; FUNCTION MIDI_DATMRG, CALDAT ; DEFINE THE KEYWORDS TO BE AVERAGED AND GET TAG NAMES OF THE STRUCT ;------------------------------------------------------------------------------- keywrd = ['CORFLX', 'MSKFLX', 'TOTFLX', 'VISAMP', 'VISPHI'] tagnam = tag_names(caldat) ; CREATE THE RESULT STRUCTURE ;------------------------------------------------------------------------------- result = caldat[0] ; LOOP OVER THE DIFFERENT KEYWORDS ;------------------------------------------------------------------------------- for i=0,n_elements(keywrd)-1 do begin ; GET THE INDICES OF THE CURRENT KEYWORD AN ITS ERROR ;----------------------------------------------------------------------- keyidx = where(tagnam eq keywrd[i]) erridx = where(tagnam eq keywrd[i]+'ERR') ; MAKE A WEIGTHED AVERAGE OF THE DATA VALUES ;----------------------------------------------------------------------- result.(keyidx) = total(caldat.(keyidx)/caldat.(erridx)^2, 2) / $ total(1./caldat.(erridx)^2, 2) ; CALCULATE THE NEW ERRORS ;----------------------------------------------------------------------- result.(erridx) = 1./sqrt(total(1./caldat.(erridx)^2, 2)) endfor ; DETERMINE THE NEW UV COORDINATES ;------------------------------------------------------------------------------- result.pbline = mean(caldat.pbline) result.pangle = mean(caldat.pangle) result.ucoord = mean(caldat.ucoord) result.vcoord = mean(caldat.vcoord) ; EDIT THE COMMENT STRING ;------------------------------------------------------------------------------- result.commen += ' MIDI_DATMRG: Average of ' + strtrim(n_elements(caldat),2) + $ ' individual data sets: ' + strjoin(caldat.identi, ', ') + '.' ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, result END ;+ ; NAME: ; MIDI_PLOT3D ; ; PURPOSE: ; Make a 3D plot of the UV plane using iPlot. ; ; CATEGORY: ; MIDI data reduction. ; ; CALLING SEQUENCE: ; midi_plot3d, caldat [, wavele] [, ranges] ; ; REQUIRED INPUTS: ; caldat a MIDI result structure containing the reduced data ; ; OPTIONAL INPUTS: ; wavele a two element array specifying the wavelength and its range ; ranges a two element array specifying the range of the U & V axes ; ; KEYWORDS: ; none ; ; OUTPUTS: ; a 3 dimensional plot of the correlated fluxes in the UV plane ; ; DESCRIPTION AND EXAMPLE: ; This procedure a 3 dimensional plot of the UV plane with the correlated ; flux in the z direction. As an input an array of structure MIDI_PRISM or ; MIDI_GRISM must be provided. The plot can then be modified and the ; visualisation changed by using further iTool commands and user inter- ; action. The baseline tracks can for example be overplotted in the UV ; plane using a command as follows, if the variable TRACK contains the U ; and V coordinates of the baseline tracks: ; midi_plot3d, caldat, [12, 0.2], [-105, 105] ; for i=0,n do iplot, track[i].ucoord, track[i].vcoord, $ ; replicate(0.0, 512), color=[200,200,200], /overplot ; ; CALLED BY: ; None. ; ; CALLING: ; None. ; ; MODIFICATION HISTORY: ; 2009-07-01 Written by Konrad R. Tristram ; PRO MIDI_PLOT3D, VISDAT, WAVELE, RANGES ; CHECK PARAMETERS ;------------------------------------------------------------------------------- if n_elements(ranges) ne 2 then ranges = [-105,105] else ranges = float(ranges) numele = n_elements(visdat) ; 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((visdat[0].wavele le wavele[0]+wavele[1]) and $ (visdat[0].wavele ge wavele[0]-wavele[1])) ; GET THE DATA TO BE PLOTTED FROM THE RESULT STRUCT ;------------------------------------------------------------------------------- z_data = reform(visdat.corflx[welmnt], n_elements(welmnt), numele) z_data = total(z_data, 1)/n_elements(welmnt) z_errs = reform(visdat.corflxerr[welmnt], n_elements(welmnt), numele) z_errs = total(z_errs, 1)/n_elements(welmnt) ; CREATE THE PLOT ;------------------------------------------------------------------------------- zrange = max(z_data + z_errs) zrange = [0, fix(zrange/10.^(fix(alog10(zrange)+10)-11)*1.1)*$ 10.^(fix(alog10(zrange)+10)-11)] iplot, [0,0], [-500,500], [0,0], xrange=-ranges, yrange=ranges, zrange=zrange, $ xtitle='U [m]', ytitle='V [m]', ztitle='correlated flux F!Dcor!N [Jy]' idplot = itgetcurrent() iplot, [-500,500], [0,0], [0,0], /overplot ; LOOP OVER THE ELEMENTS OF THE DATA AND PLOT STICKS ;------------------------------------------------------------------------------- for i=0, numele-1 do begin iplot, +visdat[i].ucoord * [1,1], +visdat[i].vcoord * [1,1], $ [0, z_data[i]], overplot=idplot, linestyle=0, $ transperancy=50.0, color=[128,128,128] iplot, -visdat[i].ucoord * [1,1], -visdat[i].vcoord * [1,1], $ [0, z_data[i]], overplot=idplot, linestyle=0, $ transperancy=50.0, color=[128,128,128] endfor ; CREATE THE SPHERE TO BE USED AS A PLOT SYMBOL ;------------------------------------------------------------------------------- sphere = obj_new('orb', color=[0,0,255], style=2, density=0.4, radius=1.0, $ shading=0) ;xobjview, sphere ; MARK THE END POINTS OF THE STICKS AND PLOT ERROR BARS ;------------------------------------------------------------------------------- iplot, [+visdat.ucoord, -visdat.ucoord], [+visdat.vcoord, -visdat.vcoord], $ [z_data, z_data], zerror=[z_errs, z_errs], overplot=idplot, $ color=[0,0,255], errorbar_color=[0,0,255], errorbar_capsize=0.02, $ linestyle=6, sym_object=sphere, thick=2 ; MARK THE BASE POINTS OF THE STICKS ;------------------------------------------------------------------------------- iplot, [+visdat.ucoord, -visdat.ucoord], [+visdat.vcoord, -visdat.vcoord], $ replicate(0.0, 2*numele), overplot=idplot, color=[0,0,255], $ linestyle=6, sym_index=7, sym_size=0.5, thick=1 END