; MODL_MKHEAD - Generate a default FITS header for the models. ; header = modl_mkhead(arrsiz, functi, object, pscale) ; ; MODL_FITS2M - Read FITS files into a model struct. ; result = modl_fits2m(patter, rotang, /silent) ; ; MODL_M2FITS - Write a model into a FITS file or into several FITS files. ; modl_m2fits, models, filnam, /cube ; ; MODL_M2VISI - Calculate the UV-plane, i.e. the complex visibilities. ; result = modl_m2visi(models, /output, /silent) ; ; MODL_M2PLOT - Plot a flux distribution and/or UV-plane to screen or PS file. ; modl_m2plot, uvplan, psfile, metafile=metafile, /flux, /fftd ; ; MODL_M2GAUS - Determine the gaussian width of a brightness distribution. ; modl_m2gaus, models [, psfile=psfile] ;+ ; NAME: ; MODL_MKHEAD ; ; PURPOSE: ; Generate a default FITS header for the models. ; ; CATEGORY: ; MIDI models. ; ; CALLING SEQUENCE: ; header = modl_mkhead(arrsiz, functi, object, pscale) ; ; REQUIRED INPUTS: ; None. ; ; OPTIONAL INPUTS: ; arrsiz 2 element array specifying the x and y size of the data ; functi string containing the name of the calling function ; object a description of the model ; pscale pixel scale of the image data ; ; KEYWORDS: ; none ; ; OUTPUTS: ; header resulting model as FITS struct ; ; DESCRIPTION AND EXAMPLE: ; This function creates a default FITS header to be used for the models. ; The function is called by: ; header = modl_mkhead([512, 512], 'MODL_GAUSSM, 'Gaussian', 1.0) ; ; CALLED BY: ; modl_gaussm ; modl_polarm ; modl_tdiskm ; modl_m2fits ; midi_2bbmdl ; midi_3bbmdl ; ; CALLING: ; Makes use of routines for FITS headers from the ASTROLIB library. ; ; MODIFICATION HISTORY: ; 2009-07-13 Written by Konrad R. Tristram ; FUNCTION MODL_MKHEAD, ARRSIZ, FUNCTI, OBJECT, PSCALE ; START GENERATING A DEFAULT HEADER FOR THE IMAGE ;------------------------------------------------------------------------------- mkhdr, header, 4, arrsiz ; ADD SEVERAL ADDITIONAL INFORMATIONAL HEADER KEYWORDS ;------------------------------------------------------------------------------- functi = strupcase(functi) sxaddpar, header, 'ORIGIN', functi, ' Created by IDL function ' + functi sxaddpar, header, 'OBJECT', object, ' Target description' sxaddpar, header, 'PSCALE', pscale, ' Pixel scale in mas per pixel' sxaddpar, header, 'WAVELE', 0.0, ' Wavelength of image in micron' sxaddpar, header, 'BUNIT ', 'Jy/Pixel', ' Units of flux: Jy per Pixel' sxaddpar, header, 'PIXSIZE', pscale/1000., ' Pixel size in arcseconds per pixel' sxaddpar, header, 'LAMBDA1', 0.0, ' Wavelength of image in m (SimVLTI)' sxaddpar, header, 'CALUNIT', 'Jy/Pix', ' Units of flux: Jy per Pixel (SimVLTI)' ; ADD AN ASTROMETRY TO THE HEADER ;------------------------------------------------------------------------------- crmatr = [[-1.0, +0.0],[+0.0, +1.0]] * pscale / 3.6E6 putast, header, crmatr, arrsiz/2, [0.0, 0.0], ['RA---TAN','DEC--TAN'] ; RESORT THE KEYWORDS AND COMMENTS IN THE HEADER ;------------------------------------------------------------------------------- header = [header[0:5], header[8:n_elements(header)-3], header[6:7], $ header[n_elements(header)-1]] ; RETURN THE HEADER ;------------------------------------------------------------------------------- return, header END ;+ ; NAME: ; MODL_FITS2M ; ; PURPOSE: ; Read FITS files into a model struct. ; ; CATEGORY: ; MIDI models. ; ; CALLING SEQUENCE: ; result = modl_fits2m(patter [, rotang], [/silent]) ; ; REQUIRED INPUTS: ; patter string containg a pattern for selecting the FITS files ; ; OPTIONAL INPUTS: ; rotang rotate the input images by this angle counterclockwise ; ; KEYWORDS: ; silent suppress output of system messages (e.g. from readfits) ; ; OUTPUTS: ; result resulting model as structure ; ; DESCRIPTION AND EXAMPLE: ; This programs takes a list of FITS files matching the pattern given ; in the input variable PATTER and reads in the data into a model for ; further analysis. The input FITS files have to all have the same ; dimensions for the image planes. Image cubes are also allowed (such ; as used e.g. by SimVLTI. The input images can be rotated by speci ; fying the rotation angle ROTANG (in degrees counterclockwise). The ; wavelength information is extracted from the header either from the ; header keyword WAVELE or assuming the SimVLTI convention: the wave- ; lengths are given sequentially for every layer in the FITS cube by ; the header keywords 'LAMBDA*'. In the latter case, the wavelengths ; are expected to be in meters, otherwise they are expected to be in ; microns. Assuming a set of fits files in the current directory called ; model08.fits, model09.fits, etc., the simplest way to run the pro- ; grammme is by calling: ; 'model = fits2model('model*.fits') ; ; CALLED BY: ; none ; ; CALLING: ; Makes use of routines for FITS headers from the ASTROLIB library. ; ; MODIFICATION HISTORY: ; 2005-10-08 Written by Konrad R. Tristram ; 2006-01-12 K. Tristram: Added support for header info on pixelscale. ; 2009-07-13 K. Tristram: Major changes including function renamed. ; 2011-02-21 K. Tristram: Fixed small inconsistency creating the struct. ; FUNCTION MODL_FITS2M, PATTER, ROTANG, SILENT=SILENT ; INITIALISATION OF PARAMETERS ;------------------------------------------------------------------------------- pscstr = 'PIXSIZE*' pscfac = 0.000 waves = 0 pxsize = 0.001 if (n_elements(rotang) eq 0) then rotang=0.0 ; RETURN TO MAIN PROGRAM LEVEL IF ERROR OCCURS ;------------------------------------------------------------------------------- ;on_error, 1 ; MAKE BACKUP OF !QUIET STATUS AND SUPPRESS SYSTEM MESSAGES IF KEYWORD IS SET ;------------------------------------------------------------------------------- quiet_status = !quiet if keyword_set(silent) then !quiet = 1 ; GET LIST OF FITS FILES MATCHING THE PATTERN ;------------------------------------------------------------------------------- if strmid(patter, strlen(patter)-5, 5) ne '.fits' then patter += '.fits' filelist = file_search(patter) ; PRINT MESSAGE AND EXIT IF NO MATCHING FITS FILES WERE FOUND ;------------------------------------------------------------------------------- if filelist[0] eq '' then begin print, 'No files found matching the file name pattern!' return, -1 endif ; START LOOPING THROUGH THE FILE LIST ;------------------------------------------------------------------------------- for i=0,n_elements(filelist)-1 do begin ; READ FITS FILE INCLUDING HEADER ;----------------------------------------------------------------------- imcube = readfits(filelist[i], header) ; GET SIZE OF THE ARRAY ;----------------------------------------------------------------------- arrsiz = size(imcube) numele = n_elements(imcube[1,1,*]) ; CHECK IF CURRENT ARRAY SIZES MATCH TO PREVIOUSLY READ DATA ;----------------------------------------------------------------------- if (n_elements(result) ne 0) then begin if (arrsiz[1] ne result[0].arrsiz[0]) or $ (arrsiz[2] ne result[0].arrsiz[1]) then begin ; PRINT A MESSAGE ;------------------------------------------------------- print, 'Dimensions of the input FITS files ', $ "don't match! Exiting with current data only!" ; RESET QUIET SYSTEM VARIABLE AND RETURN THE RESULT ;------------------------------------------------------- !quiet = quiet_status return, result endif endif ; FIND OUT IN WHICH HEADER KEYWORDS THE WAVELENGTHS ARE GIVEN ;----------------------------------------------------------------------- idxvar = where(strmatch(header, 'LAMBDA*', /fold)) if (n_elements(idxvar) ne numele) then begin idxvar = where(strmatch(header, 'WAVELE*', /fold)) tmpvar = 1.0 ; wavelengths should be in micron endif else tmpvar = 1.0E+06 ; wavelengths should be in meters ; IF THERE ARE AS MANY WAVELENGTHS AS ARE IMAGE PLANES, THEN GET THEM ;----------------------------------------------------------------------- if (n_elements(idxvar) eq numele) and (idxvar[0] ne -1) then begin wavele = float(strmid(header[idxvar], 11, 20)) * tmpvar ; OTHERWISE PRINT A MESSAGE AND SET A DEFAULT VALUE ;----------------------------------------------------------------------- endif else begin print, 'Something is wrong with wavelength information!' wavele = findgen(numele) endelse ; FOR THE PIXEL SCALE, CHECK IF THE FITS HEADER CONTAINS AN ASTROMETRY ;----------------------------------------------------------------------- extast, header, astrom, numpar if numpar ge 1 then begin ; GET THE PIXEL SCALE FROM THE ASTROMETRY ;--------------------------------------------------------------- getrot, astrom, tmpvar, pscale ; CHECK IF PIXEL SCALE IS SYMMETRIC ;--------------------------------------------------------------- if abs(pscale[0]) eq abs(pscale[1]) then begin pscale = abs(pscale[0]) * 1000 * 3600 endif else begin print, ' Can not work with different pixel scales' + $ ' in x and y direction.' pscale = -1.0 endelse ; OTHERWISE TRY TO READ THE THE PIXEL SCALE FROM THE KEYWORD 'PIXSIZE' ;----------------------------------------------------------------------- endif else begin idxvar = where(strmatch(header, 'PIXSIZE*', /fold)) if idxvar[0] ge 0 then begin pscale = abs(float(strmid(header[idxvar], 11, 20))) pscale *= 1000. ; pixel scale should be in arcsec / pix endif else begin print, 'No pixel scale defined! Assuming 1 mas / pixel' pscale = 1.0 endelse endelse ; GET OBJECT KEYWORD FROM THE FITS HEADER ;----------------------------------------------------------------------- object = sxpar(header, 'OBJECT') if size(object, /type) ne 7 then object = 'Model read from ' + $ filelist[i] ; LOOP THROUGH PLANES OF CURRENT FITS CUBE ;----------------------------------------------------------------------- for j=0,numele-1 do begin ; ROTATE IMAGE BY SPECIFIED ANGLE ;--------------------------------------------------------------- imgdat = rot(imcube[*,*,j], -rotang, /interp) ; DETERMINE SIZE OF ARRAY AND CALCULATE NEW SIZE ;--------------------------------------------------------------- newsiz = min([arrsiz[1], arrsiz[2]]) / $ (abs(cos(rotang/180.*!pi)) + $ abs(sin(rotang/180.*!pi))) newsiz = fix(newsiz+0.5) ; CROP FLUX ARRAY TO NEW SIZE ;--------------------------------------------------------------- imgdat = imgdat[arrsiz[1]/2-newsiz/2:arrsiz[1]/2+newsiz/2-1, $ arrsiz[2]/2-newsiz/2:arrsiz[2]/2+newsiz/2-1] ; UPDATE THE FITS HEADER KEYWORDS ;--------------------------------------------------------------- check_fits, imgdat, header, /update sxaddpar, header, 'WAVELE', wavele[j], ' Wavelength of ' + $ 'image in micron', after='PSCALE' sxaddpar, header, 'LAMBDA1', wavele[j]*1.0E-6, $ ' Wavelength of image in m (SimVLTI)', after='PIXSIZE' ;#################### TESTING!!! putast, header, [[1,0], [0,1]]*pscale[0]/3600./180./1000, $ arrsiz[1:2]/2, [0,0], ['RA---TAN','DEC--TAN'] ;#################### TESTING!!! ; CREATE FIRST ELEMENT OF RESULTING STRUCT ;--------------------------------------------------------------- newsiz *= [1,1] if (n_elements(result) eq 0) then begin result = {object:object, wavele:wavele[j], $ arrsiz:newsiz, pscale:pscale[0], $ header:ptr_new(header), imgdat:imgdat} ; OR APPEND THIS PLANE TO THE EXISTING STRUCT ;--------------------------------------------------------------- endif else begin result = [result, {object:object, wavele:wavele[j], $ arrsiz:newsiz, pscale:pscale[0], $ header:ptr_new(header), imgdat:imgdat}] endelse endfor endfor ; RESET QUIET SYSTEM VARIABLE TO ITS ORIGINAL STATE AND RETURN THE RESULT ;------------------------------------------------------------------------------- !quiet = quiet_status return, result END ;+ ; NAME: ; MODL_M2FITS ; ; PURPOSE: ; Write a model into a FITS file or into several FITS files. ; ; CATEGORY: ; MIDI models. ; ; CALLING SEQUENCE: ; modl_m2fits, models [, filnam] [, /cube] ; ; REQUIRED INPUTS: ; models model struct to be written into the FITS file. ; ; OPTIONAL INPUTS: ; filnam path and name of the file where the data should ; be written to, e.g. filename = '/tmp/test.fits' ; ; KEYWORDS: ; cube write out the models into a single FITS cube ; ; OUTPUTS: ; - a FITS file is written into the designated file. ; ; DESCRIPTION AND EXAMPLE: ; This function writes a model struct into a single FITS file or into ; several FITS files. The input model must be a struct of the form ; {wavele:0.0, pscale:0.0, imgdat:fltarr(n,n)}, such as created by e.g. ; 'models = modl_gaussm()'. When the keyword CUBE is set, then a single ; FITS cube is created (assuming the same pixel scale in all images), ; otherwise a set of FITS images is produced. With CUBE set, the ; resulting file is compatible with SimVLTI. The programme is simply ; called by 'model2fits, model'. ; ; CALLED BY: ; none ; ; CALLING: ; Makes use of routines for FITS headers from the ASTROLIB library. ; modl_mkhead ; ; MODIFICATION HISTORY: ; 2005-09-05 Written by Konrad R. Tristram ; 2009-07-13 K. Tristram: Major changes including function renamed. ; PRO MODL_M2FITS, MODELS, FILNAM, CUBE=CUBE ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- if (n_elements(filnam) eq 0) then filnam = 'geomodel' numele = n_elements(models) ; CHECK IF A HEADER IS PRESENT IN THE STRUCT ;------------------------------------------------------------------------------- if total(tag_names(models[0]) eq 'HEADER') then begin check_fits, models[0].imgdat, models[0].header, errmsg=errmsg endif else errmsg = 'No FITS header present in the struct.' ; CREATE A DEFAULT HEADER IF NONE IS PRESENT IN THE STRUCTURE ;------------------------------------------------------------------------------- if errmsg ne '' then begin header = modl_mkhead(arrsiz, 'MODL_M2FITS', '', models[0].pscale) endif else header = models[0].header ; CREATE A SINGLE FITS FILE WITH THE DATA IN A FITS CUBE ;------------------------------------------------------------------------------- if keyword_set(cube) then begin ; UPDATE THE FITS HEADER KEYWORDS AND REMOVE THE KEYWORD 'WAVELE' ;----------------------------------------------------------------------- check_fits, models.imgdat, header, /update sxdelpar, header, 'WAVELE' ; CHECK IF "OBJECT" TAG NAME IS DEFINED IN THE MODEL STRUCTURE ;----------------------------------------------------------------------- if total(tag_names(models) eq 'OBJECT') then begin object = models.object ; PRINT WARNING IF NOT ALL DESCRIPTIONS ARE IDENTICAL ;--------------------------------------------------------------- if ~ array_equal(object, shift(object, 1)) then begin print, 'Warning: different model descriptions! ' + $ 'Only using the first one in the FITS file.' endif ; SAVE THE FIRST DESCRIPTION INTO THE FITS HEADER ;--------------------------------------------------------------- sxaddpar, header, 'OBJECT', object[0], before='PSCALE' ; OTHERWISE SAVE ALTERNATIVE OBJECT DESCRIPTION INTO THE HEADER ;----------------------------------------------------------------------- endif else sxaddpar, header, 'OBJECT', 'Model written by MODL_M2FITS', $ before='PSCALE' ; LOOP OVER THE WAVELENGTHS AND ADD WAVELENGTH KEYWORDS TO THE HEADER ;----------------------------------------------------------------------- for i=0,numele-1 do begin sxaddpar, header, 'LAMBDA' + strcompress(i+1, /remove_all), $ models[i].wavele*1.E-6, before='CALUNIT' endfor ; WRITE THE FITS CUBE ;----------------------------------------------------------------------- writefits, filnam+'.fits', models.imgdat, header ; OTHERWISE CREATE A FITS FILE FOR EVERY ELEMENT IN MODELS ;------------------------------------------------------------------------------- endif else begin for i=0, numele-1 do begin ; CHECK THE FITS HEADER KEYWORDS ;--------------------------------------------------------------- header = models[i].header check_fits, models[i].imgdat, header, /update ; CHECK IF "OBJECT" TAG NAME IS DEFINED IN THE MODEL STRUCTURE ;--------------------------------------------------------------- if total(tag_names(models) eq 'OBJECT') then begin ; SAVE THE MODEL DESCRIPTION INTO THE FITS HEADER ;------------------------------------------------------- sxaddpar, header, 'OBJECT', models[i].object, $ before='PSCALE' ; OTHERWISE SAVE ALTERNATIVE OBJECT DESCRIPTION INTO THE HEADER ;--------------------------------------------------------------- endif else sxaddpar, header, 'OBJECT', 'Model written by ' + $ 'MODL_M2FITS', before='PSCALE' ; WRITE THE FITS FILE ;--------------------------------------------------------------- ;wavstr = string(models[i].wavele, format='(I02, ".fits")') wavstr = string(i, format='(I02, ".fits")') writefits, filnam+wavstr, models[i].imgdat, header endfor endelse END ;+ ; NAME: ; MODL_M2VISI ; ; PURPOSE: ; Calculate the UV-plane, i.e. the complex visibilities for a model. ; ; CATEGORY: ; MIDI models. ; ; CALLING SEQUENCE: ; result = modl_m2visi(models [, /nophas] [, /output] [, /silent]) ; ; REQUIRED INPUTS: ; models model struct for which the UV-plane is to be calculated. ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORDS: ; output make plots of the flux distribution and UV-plane ; silent supress output of processor time ; ; OUTPUTS: ; result resulting UV-plane as a structure ; ; DESCRIPTION AND EXAMPLE: ; This program calculates the UV-plane for a model by simply carrying ; out a fast fourier transform. The input model must be a struct of the ; form {wavele:0.0, pscale:0.0, imgdat:fltarr(n,n)}, such as created by ; e.g. 'models = modl_gaussm()'. The result contains the complex result ; of the Fourier Transform. The programme can the be simply called by: ; result = model2vis(model) ; ; CALLED BY: ; midi_2bbmdl ; midi_3bbmdl ; ; CALLING: ; Makes use of routines for FITS headers from the ASTROLIB library. ; modl_mkhead ; modl_m2plot ; ; MODIFICATION HISTORY: ; 2005-09-08 Written by Konrad R. Tristram ; 2005-09-09 K. Tristram: Added calculation of computation time. ; Moved plot part to back of program. ; 2005-10-06 K. Tristram: Removed errors in documentation. ; 2006-01-13 K. Tristram: Added support for pixelscale. ; 2006-06-26 K. Tristram: Made RA axis count in correct direction in plot. ; 2006-12-06 K. Tristram: Moved plotting code to extra plotting routine. ; 2007-02-07 K. Tristram: Added the extraction of phases. ; 2009-07-13 K. Tristram: Function renamed and rewritten. ; FUNCTION MODL_M2VISI, MODELS, OUTPUT=OUTPUT, SILENT=SILENT ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- systim = systime(/seconds) ; GENERATE STRUCT HOLDING THE RESULT ;------------------------------------------------------------------------------- arrsiz = size(models[0].imgdat) result = create_struct(models[0], 'fftdat', complexarr(arrsiz[1:2])) result = replicate(result, n_elements(models)) ; COPY THE INPUT MODEL INTO THE RESULT STRUCTURE ;------------------------------------------------------------------------------- for i=0,n_tags(models)-1 do result.(i) = models.(i) ; LOOP OVER WAVELENGTH ELEMENTS IN THE MODEL ;------------------------------------------------------------------------------- for i=0, n_elements(models)-1 do begin ; GET THE SIZE OF THE IMAGE ;----------------------------------------------------------------------- sz_img = size(models[i].imgdat) ; DO THE FOURIER TRANSFORM AND SHIFT IT ;----------------------------------------------------------------------- fftdat = fft(shift(models[i].imgdat, -sz_img[1]/2, -sz_img[2]/2)) fftdat = shift(fftdat, sz_img[1]/2, sz_img[2]/2) ; NORMALISE THE RESULT AND COPY IT INTO THE RESULT STRUCTURE ;----------------------------------------------------------------------- factor = total(models[i].imgdat) / max(abs(fftdat)) result[i].fftdat = fftdat * factor endfor ; PRINT TIME PROGRAMME TOOK TO RUN ;------------------------------------------------------------------------------- if ~ keyword_set(silent) then begin print, 'Computation time: ', systime(/seconds) - systim, ' sec' endif ; SHOW RESULT IF DESIRED ;------------------------------------------------------------------------------- if keyword_set(output) then modl_m2plot, result ; SUPPRESS OUTPUT OF MATH ERRORS ;------------------------------------------------------------------------------- dump = check_math() ; RETURN RESULT ;------------------------------------------------------------------------------- return, result END ;+ ; NAME: ; MODL_M2PLOT ; ; PURPOSE: ; Plot a model flux distribution and/or UV-plane to screen or PS file. ; ; CATEGORY: ; MIDI models. ; ; ; CALLING SEQUENCE: ; modl_m2plot, uvplan [, psfile] [, metafile=metafile] [, minval=minval] ; [, /flux] [, /fftd] [, delay=delay], [, factor=factor] ; [, maxpbl=maxpbl] [, overpl=overpl] ; ; REQUIRED INPUTS: ; uvplan struct containing the data to be plotted ; ; OPTIONAL INPUTS: ; psfile switches to PS and sets filename base for the PS file ; ; KEYWORDS: ; metafile switches to WMF output and sets the filename base ; minval minimum value for surface plot of the flux ; flux switch for plot of the flux distribution ; fftd switch for plot of the UV-plane ; delay delay in seconds between creating the plots ; factor magnification factor for the plot ; maxpbl maximum projected baseline length to plot in a UV plot ; overpl do not open a new window but reuse the current one ; ; OUTPUTS: ; - plots of the flux or the UV plane as desired ; ; DESCRIPTION AND EXAMPLE: ; This program plots out either the flux distribution or the Fourier ; transform of the flux distribution (the UV-plane) or both at once. ; The input for this has to be a struct containing the corresponding ; data, e.g. as created by 'model=gaussmodel()' or by ; 'uvplan = model2vis(model)'. The plot can displayed on the screen, or ; a PS file can be created when the parameter PSFILE is set. ; ; CALLED BY: ; modl_gaussm ; modl_polarm ; modl_tdiskm ; modl_m2visi ; midi_extvis ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2006-09-08 Written by Konrad R. Tristram (Adapted from MODEL2VIS) ; 2006-12-18 K. Tristram: Added support for WMF output. ; 2007-06-10 K. Tristram: Set tick marks to the outside in PS plot. ; 2009-07-13 K. Tristram: Function renamed and several keywords added. ; 2012-11-01 K. Tristram: Adjust PS font size for !p.font=0. ; PRO MODL_M2PLOT, UVPLAN, PSFILE, METAFILE=METAFILE, MINVAL=MINVAL, FLUX=FLUX, $ FFTD=FFTD, DELAY=DELAY, FACTOR=FACTOR, MAXPBL=MAXPBL, $ OVERPL=OVERPL ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- if (n_elements(minval) eq 0) then minval = 1.0e-20 if (n_elements(factor) eq 0) then factor = 1.0 if (n_elements(flux) eq 0) then flux=byte(1) else flux=byte(flux) if (n_elements(fftd) eq 0) then fftd=byte(1) else flux=byte(flux) if not total(tag_names(uvplan) eq 'IMGDAT' or $ tag_names(uvplan) eq 'WAVELE') then begin print, 'No flux of fftd information present in structure!' return endif if not total(size(delay, /type) eq [1,2,3,4,5,7]) then delay = 0.5 $ else delay = delay[0] ; BACKUP PLOT VARIABLE ;------------------------------------------------------------------------------- P = !P ; GET MIN AND MAX OF THE FLUX FOR OVERALL SCALING ;------------------------------------------------------------------------------- if total(tag_names(uvplan) eq 'IMGDAT') then begin flxmax = alog(max(uvplan.imgdat > minval)) flxmin = alog(min(uvplan.imgdat > minval)) sz_img = size(uvplan.imgdat) if (fftd and total(tag_names(uvplan) eq 'FFTDAT')) eq 0 then flux=1B endif else flux=byte(0) if total(tag_names(uvplan) eq 'FFTDAT') then begin fftabs = abs(uvplan.fftdat) fftmax = alog(max(fftabs > minval)) fftmin = alog(min(fftabs > minval)) sz_img = size(fftabs) if flux eq 0 then fftd=1B endif else fftd=byte(0) ; OPEN DEVICE, CHANGE DEVICE SETTINGS AND SET COORDINATES FOR PLOTTING ;------------------------------------------------------------------------------- if keyword_set(psfile) then begin d_name = !d.name set_plot, 'ps' imgsiz = sz_img[1] pltsiz = 8000 ps_xs = 11 ps_xo = 5 !p.charsize = 1.0 !p.ticklen = -0.02 if (flux + fftd) eq 2 then begin posit1 = [2200,11100] posit2 = [2200, 1100] posit3 = [ 0,10000] ps_ys = 20 ps_yo = 6 endif else begin posit1 = [2200, 1500] posit2 = [2200, 1500] posit3 = [ 0,10000] ps_ys = 11 ps_yo = 11 endelse ;For publications take these formats by uncommenting the next lines ps_xs=8.0 & ps_ys=8.0 & ps_xo=6.5 & ps_yo=10.5 & pltsiz = 6000 posit1=[1500, 1100] & posit3 = [ 0, 7500] if !p.font eq 0 then !p.charsize=0.65 else !p.charsize=0.75 endif else if keyword_set(metafile) then begin ;For POWERPOINT then have to rescale plot to 25% d_name = !d.name set_plot, 'metafile' !p.font=0 !p.charsize=0.80 imgsiz = sz_img[1] pltsiz = 1700 ps_xs = 40 if (flux + fftd) eq 2 then begin posit1 = [220,2180] posit2 = [220, 180] posit3 = [ 0,2000] ps_ys = 80 endif else begin posit1 = [220, 180] posit2 = [220, 180] posit3 = [100,2000] ps_ys = 40 endelse endif else begin if !version.os_family eq 'Windows' then begin set_plot, 'WIN' endif else set_plot, 'X' ;device, decomposed=0 pltsiz = 450 * factor imgsiz = 450 * factor posit2 = [55, 35] posit3 = [60, pltsiz + 45] !p.charsize=1.05 if ~ keyword_set(overpl) then begin if (flux + fftd) eq 2 then begin window, /FREE, xsize=pltsiz+70, ysize=2*pltsiz+110 posit1 = [55, pltsiz + 85] endif else begin window, /FREE, xsize=pltsiz+70, ysize=pltsiz+60 posit1 = [55, 35] endelse endif else erase endelse ; LOAD COLOR TABLE "HEAT" ;------------------------------------------------------------------------------- loadct, 3, /silent ; LOOP OVER WAVELENGTH ELEMENTS ;------------------------------------------------------------------------------- for i=0, n_elements(uvplan)-1 do begin ; OPEN PLOT DEVICE ;----------------------------------------------------------------------- if keyword_set(psfile) then device, file = psfile + $ string(uvplan[i].wavele*10, format='(I03)')+'.ps', $ xsize=ps_xs, ysize=ps_ys, xoffset=ps_xo, yoffset=ps_yo, $ bits_per_pixel=8, /portrait, /color $ else if keyword_set(metafile) then device, file = metafile + $ string(uvplan[i].wavele*10, format='(I03)')+'.emf', $ set_font='Arial*76', xsize=ps_xs, ysize=ps_ys ; MAKE SURFACE PLOT OF FLUX DISTRIBUTION ;----------------------------------------------------------------------- if flux then begin ; RESIZE AND RESCALE THE ARRAY: 1. LOGARITHMIC ;--------------------------------------------------------------- pltarr = congrid(alog(uvplan[i].imgdat > minval), imgsiz, imgsiz) pltarr = (pltarr-flxmin)/(flxmax-flxmin) * 255 ; RESIZE AND RESCALE THE ARRAY: 2. SQUARE ROOT ;--------------------------------------------------------------- ; pltarr = congrid(sqrt(uvplan[i].imgdat), imgsiz, imgsiz) ; pltarr = (pltarr / 0.04 * 235. < 235.) + 20. ; MAKE PLOT OF THE ARRAY ;--------------------------------------------------------------- tv, pltarr, posit1[0], posit1[1], $ xsize=pltsiz, ysize=pltsiz, /DEVICE ; PLOT AXES AND LINES ON TOP OF THE ARRAY ;--------------------------------------------------------------- ranges = sz_img[1] * uvplan[i].pscale / 2. plot,[0,0], xrange=[ ranges, -ranges], $ yrange=[-ranges, ranges], $ xstyle=1, ystyle=1, position=[posit1, posit1+pltsiz], $ xtitle='RA in mas', ytitle='DEC in mas', $ title='Flux distribution', $ /DEVICE, /NODATA, /NOERASE oplot, [0,0], [-1000, 1000] oplot, [-1000, 1000], [0,0] endif ; MAKE SURFACE PLOT OF FOURIER TRANSFORM ;----------------------------------------------------------------------- if fftd then begin pltarr = fftabs[*,*,i] if n_elements(maxpbl) ge 1 then begin pixran = maxpbl[0] * !pi/(uvplan[i].wavele * 1.0e-6) / $ (3600.*180.*1000.) * sz_img[1]*uvplan[i].pscale pixran = pixran < sz_img[1]/2-1 pltarr = pltarr[sz_img[1]/2-pixran:sz_img[1]/2+pixran, $ sz_img[2]/2-pixran:sz_img[2]/2+pixran] endif tmpsiz = size(pltarr) ; RESIZE AND RESCALE THE ARRAY ;--------------------------------------------------------------- pltarr = congrid(alog(pltarr > minval), imgsiz, imgsiz) pltarr = (pltarr-fftmin)/(fftmax-fftmin) * 255 ; MAKE PLOT OF THE ARRAY ;--------------------------------------------------------------- tv, pltarr, posit2[0], posit2[1], $ xsize=pltsiz, ysize=pltsiz, /DEVICE ; PLOT AXES AND LINES ON TOP OF THE ARRAY ;--------------------------------------------------------------- ranges = (uvplan[i].wavele * 1.0e-6) / !pi * $ (3600. * 180. * 1000.) / uvplan[i].pscale * $ tmpsiz[1] / sz_img[1] if not (keyword_set(psfile) or keyword_set(metafile)) then begin tv, make_array(pltsiz+70, 28, value=1.0), 0, 0 tv, make_array(53, pltsiz+50, value=1.0), 0, 0 endif plot,[0,0], xrange=[ ranges/2., -ranges/2.], $ yrange=[-ranges/2., ranges/2.], $ xstyle=1, ystyle=1, position=[posit2, posit2+pltsiz], $ xtitle='U in m', ytitle='V in m', $ title='Fourier transform (UV-Plane)', $ /DEVICE, /NODATA, /NOERASE oplot, [0,0], [-100000, 100000] oplot, [-100000, 100000], [0,0] endif ; WRITE OUT CURRENT WAVELENGTH ;----------------------------------------------------------------------- if not (keyword_set(psfile) or keyword_set(metafile)) then begin tv, make_array(100, 20, value=1.0), 50, pltsiz+40 endif xyouts, posit3[0], posit3[1], $ string(uvplan[i].wavele,format='(F5.1)')+' !4l!xm', $ charsize=1.3*!p.charsize, /DEVICE if keyword_set(psfile) then device, /close ; WAIT FOR USER INPUT OR SOME PREDEFINED PERIOD BEFORE THE NEXT FRAME ;----------------------------------------------------------------------- if size(delay, /type) eq 7 then begin read, 'Press ENTER to continue, "f" to finish!', delay on_ioerror, nonumbr if delay eq 'f' then delay = 0.0 else $ if float(delay) ne 0.0 then delay = float(delay) nonumbr: endif else wait, delay endfor ; CHANGE BACK DEVICE SETTINGS AND LOAD DEFAULT COLOUR MAP ;------------------------------------------------------------------------------- if keyword_set(psfile) or keyword_set(metafile) then begin device, /close set_plot, d_name endif loadct, 0, /silent ; RESTORE PLOT VARIABLE ;------------------------------------------------------------------------------- !P = P END ;+ ; NAME: ; MODL_M2GAUS ; ; PURPOSE: ; Determine the gaussian width of a projected brightness distribution. ; ; CATEGORY: ; MIDI Analysis. ; ; CALLING SEQUENCE: ; modl_m2gaus, models [, psfile=psfile] ; ; REQUIRED INPUTS: ; models a model struct containing a brightness distribution ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORDS: ; psfile name for PS file and trigger for PS output ; ; OUTPUTS: ; - plot of fitted gauss widths on screen or into PS file ; ; DESCRIPTION AND EXAMPLE: ; This procedure uses a model brightness distribution as input and ; collapses the distribution along a certain direction. The one ; dimensional distribution is then fitted by a Gaussian distribution. ; This procedure is repeated for a number of position angles. The FWHM ; of this Gaussian is plotted in a polar plot in the top left, a plot of ; the FWHM depending on the position angle is shown in the lower left. ; The original flux distribution is depicted in the right plot. The input ; model has to be a struct containing at least a field containg the flux. ; The simplest way to call the procedure is: 'modl_m2gaus, model'. ; ; MODIFICATION HISTORY: ; 2006-05-01 Written by Konrad R. Tristram ; 2006-06-26 K. Tristram: Made RA axis count in correct direction in plot. ; 2009-07-13 K. Tristram: Function renamed plus several minor changes. ; PRO MODL_M2GAUS, MODELS, PSFILE=PSFILE ; GET SIZE OF THE INPUT ARRAY ;------------------------------------------------------------------------------- x_size = n_elements(models.imgdat(*,0)) / 2. y_size = n_elements(models.imgdat(0,*)) / 2. ; SET SAMPLING IN ANGULAR DIRECTION ;------------------------------------------------------------------------------- angles = findgen(121)*3.0 ; SET PLOT AND TEXT POSITIONS ;------------------------------------------------------------------------------- plotpos1 = [ 30,240,390,590] plotpos2 = [ 30, 30,390,210] plotpos3 = [430, 30,990,590] textposi = [[450, 570], [560, 560]] textoffs = [[ 0, 0], [-20, -20]] ; GENERATE ARRAY CONTAING FWHM OF GAUSSIAN FITS ;------------------------------------------------------------------------------- fwhm = make_array(n_elements(angles), /float, value=0.0) ; LOOP ON ANGLES ;------------------------------------------------------------------------------- for i=0,n_elements(angles)-1 do begin ; COLLAPSE FLUX DISTRIBUTION ;----------------------------------------------------------------------- colaps = total(rot(models[0].imgdat, -angles[i]), 1) ; DO GAUSSFIT AND SAVE RESULT IN VARIABLE ;----------------------------------------------------------------------- result = gaussfit(findgen(n_elements(colaps)), colaps, $ params, nterms=3) fwhm(i) = params(2)* 2 * sqrt(2*alog(2)) endfor ; SET PLOT DEVICE AND MAKE DEVICE SPECIFIC SETTINGS ;------------------------------------------------------------------------------- if keyword_set(psfile) then begin set_plot, 'ps' device, file=psfile, xsize=21, ysize=13, bits_per_pixel=8, /color plotpos1 = plotpos1 * 20.0 + 400 plotpos2 = plotpos2 * 20.0 + 400 plotpos3 = plotpos3 * 20.0 + 400 textposi = textposi * 20.0 + 400 textoffs = textoffs * 20.0 color1 = 240 color2 = 200 charsz = 0.6 endif else begin window, /FREE, xsize=1000, ysize=600 color1 = 240 color2 = 200 charsz = 1.0 endelse ; GENERATE POLAR PLOT OF FWHM, PLOT CROSS IN THE MIDDLE ;------------------------------------------------------------------------------- loadct, 39, /silent x = fwhm * sin(angles/180*!pi) / 2.0 y = fwhm * cos(angles/180*!pi) / 2.0 plot, [0, 0], [0, 0], xrange=[0.6*max(fwhm), -0.6*max(fwhm)], $ yrange=[-0.6*max(fwhm), 0.6*max(fwhm)], $ position=plotpos1, charsize=charsz, /DEVICE, /NODATA oplot, -x, y, color=color1 plots, 0, 0, psym=1 ; PLOT FWHM VERSUS ANGLE ;------------------------------------------------------------------------------- plot, [0, 0], [0, 0], xrange=[min(angles), max(angles)], xstyle=1, $ yrange=[0, 1.2*max(fwhm)], position=plotpos2, charsize=charsz, $ /DEVICE, /NODATA, /NOERASE oplot, angles, fwhm, color=color1 ; MAKE GRAYSCALE PLOT OF ORIGINAL FLUX DISTRIBUTION ;------------------------------------------------------------------------------- pltarr = models.imgdat flxmax = max(pltarr) flxmin = min(pltarr) loadct, 0, /silent if keyword_set(psfile) then begin pltarr = 255 - (pltarr-flxmin)/(flxmax-flxmin) * 255 tv, pltarr, plotpos3[0], plotpos3[1], xsize=plotpos3[2]-plotpos3[0], $ ysize=plotpos3[3]-plotpos3[1], /DEVICE endif else begin pltarr = congrid(pltarr, 560, 560) pltarr = (pltarr-flxmin)/(flxmax-flxmin) * 255 tv, pltarr, plotpos3[0], plotpos3[1] endelse ; GENERATE AXIS AROUND GRAYSCALE PLOT AND PLOT CROSS IN THE MIDDLE ;------------------------------------------------------------------------------- loadct, 39, /silent plot,[0,0], xrange=[x_size, -x_size], yrange=[-y_size, y_size], xstyle=1, $ ystyle=1, position=plotpos3, charsize=charsz, /DEVICE, /NODATA, /NOERASE plots, 0, 0, psym=1, color=color2 ; GENERATE ARRAY CONTAINING KNOWN PARAMETERS AND GET NAMES IN MODEL STRUCTURE ;------------------------------------------------------------------------------- tagstr = ['FWHMAX', 'IANGLE', 'PANGLE', 'EXPONE', 'PORDER'] tagnam = tag_names(models) ; LOOP ON KNOWN PARAMETERS ;------------------------------------------------------------------------------- for i=0,n_elements(tagstr)-1 do begin ; DETERMINE IF STRUCT CONTAINS KNOWN PARAMETER ;----------------------------------------------------------------------- tagpos = where(strcmp(tagnam, tagstr[i])) if (tagpos ge 0) then begin ; PRINT OUT PARAMETER AND ITS VALUE ;--------------------------------------------------------------- xyouts, textposi[0,0], textposi[0,1], tagstr[i] + ': ', $ color=color1, charsize=1.2*charsz, /DEVICE xyouts, textposi[1,0], textposi[1,1], string(models.(tagpos), $ format='(F7.2)'), align=1.0, color=color1, $ charsize=1.2*charsz, /DEVICE ; GO TO NEW LINE FOR NEXT TEXT OUTPUT ;--------------------------------------------------------------- textposi = textposi + textoffs endif endfor ; CHANGE BACK DEVICE SETTINGS AND LOAD DEFAULT COLOUR MAP ;------------------------------------------------------------------------------- if keyword_set(psfile) then begin device, /close set_plot, 'X' loadct, 0, /silent endif END