; MIDI_2BBPRE - Example routine: Create the variables needed for a fit. ; transf = midi_2bbpre() ; ; MIDI_2BBMDL - Create a 2 component 2D Gaussian black body brightness profile. ; gmodel = midi_2bbmdl(arrsiz, wavele, params, tausil=tausil) ; ; MIDI_2BBDAT - Calculate and return the "visibility data" for the model. ; mdlvis = midi_2bbdat(transf, params, /clumps, /reclmp) ; ; MIDI_2BBFIT - Run the fit for the 2 component Gaussian black body model. ; midi_2bbfit, transf, result ; ; MIDI_2BBCMP - Compare a 2 component Gaussian black body model to the data. ; midi_2bbcmp, transf, params, wavele, psfile='filename' ;+ ; NAME: ; MIDI_2BBPRE ; ; PURPOSE: ; Example routine: Create the variables needed for a fit. ; ; CATEGORY: ; MIDI two component black body fit. ; ; CALLING SEQUENCE: ; transf = midi_2bbpre() ; ; REQUIRED INPUTS: ; none ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; none ; ; OUTPUTS: ; transf transfer variable containing information on data and model ; ; DESCRIPTION AND EXAMPLE: ; This is an example for a function to generate the necessary struc- ; ture needed to start a fit of MIDI data with MIDI_2BBFIT. First a ; few general parameters of the fit are defined and the silicate ab- ; sorption profile is loaded. Then the MIDI data is loaded from hard ; disk and reorganised so that it can be used for the modelling. All ; the necessary options for the fitting have to be set using this ; function. The return value of the function is a structure which ; contains the information needed for the generation of the model ; and the fitting procedure. The calling sequence is: ; transf = midi_2bbpre() ; ; CALLED BY: ; none ; ; CALLING: ; midi_prepht ; ; MODIFICATION HISTORY: ; 2007-02-07 Adapted from earlier test_5.pro by Konrad R. Tristram ; 2009-07-22 K. Tristram: Renamed and adapted to the new fitting routine. ; FUNCTION MIDI_2BBPRE ; SET THE ARRAY SIZE AND PIXEL SCALE FOR THE MODEL ; ------------------------------------------------------------------------------ arrsiz = 512 pscale = 2 ; SELECT WHICH DATA TO USE FOR THE FITTING ; ------------------------------------------------------------------------------ ; corflx, visamp, visphi select = [ 1, 0, 0] ; fit correlated fluxes and phases ; SET THE WAVELENGTH RANGE AND SPACING FOR THE COMPARISON ; ------------------------------------------------------------------------------ wavele = findgen(21)/4. + 8.0 ; use 21 wavelengths from 8 to 13 micron ; RESTORE THE SILICATE ABSORPTION PROFILE (OPTICAL DEPTH) ; ------------------------------------------------------------------------------ restore, 'tausil.sav' ; RESTORE THE DATA ; ------------------------------------------------------------------------------ restore, 'mydata.sav' ; SELECT THE UV POINTS TO BE USED FOR THE DETERMINATION OF THE TOTAL FLUX ; ------------------------------------------------------------------------------ selflx = indgen(n_elements(mydata)) ; simply take all photometries ; CALCULATE THE PHOTOMETRY AND PREFIX IT AS BASELINE WITH LENGTH ZERO ; ------------------------------------------------------------------------------ datvis = midi_prepht(mydata, selflx) ; SELECT WHICH UV POINTS TO USE FOR THE FITTING ; ------------------------------------------------------------------------------ seluvp = indgen(n_elements(datvis)) ; select all visibility points measured ; DO THE SELECTION OF THE BASELINES AND THEN SORT THEM BY POSITION ANGLE ; ------------------------------------------------------------------------------ datvis = datvis[seluvp] datvis = datvis[reverse(sort(datvis.pangle))] datvis = datvis[reverse(sort(datvis.baseli))] ; SET THE STARTING PARAMETERS AND WHICH PARAMETERS TO FIT ; ------------------------------------------------------------------------------ ; x-off1 y-off1 fwhm1 iangl1 pangl1 sildp1 temp1 fact1 x-off2 y-off2 fwhm2 iangl2 pangl2 sildp2 temp2 fact2 params = [ 0.000, 0.000, 15.00, 60.00, 50.00, 2.000, 500.0, 1.000, 10.00, -5.00, 80.00, 30.00, 10.00, 1.000, 300.0, 0.500] ; generator ;params = [ 0.000, 0.000, 14.94, 58.36, 47.75, 1.869, 514.8, 0.843, 10.30, -4.92, 79.90, 30.04, 9.242, 0.915, 292.8, 0.539] ; Xr^2 = 0.313 params = [ 0.000, 0.000, 20.00, 80.00, 0.00, 1.500, 400.0, 1.000, 0.000, 0.000, 100.0, 20.00, 0.00, 1.500, 250.0, 1.000] ; random guess ;params = [ 0.000, 0.000, 14.84, 59.24, 47.73, 1.900, 510.5, 0.891, 12.08,-0.037, 80.06, 31.36, 10.87, 0.927, 296.4, 0.519] ; Xr^2 = 0.346 parfix = [ 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] ; RETURN THE STRUCT NEEDED FOR THE MODEL FITTING ; ------------------------------------------------------------------------------ return, {arrsiz:arrsiz, wavele:wavele, pscale:pscale, params:params, $ parfix:parfix, tausil:tausil, datvis:datvis, select:select} END ;+ ; NAME: ; MIDI_2BBMDL ; ; PURPOSE: ; Create a 2 component 2D Gaussian black body brightness distribution. ; ; CATEGORY: ; MIDI two component black body fit. ; ; CALLING SEQUENCE: ; gmodel = midi_2bbmdl([arrsiz] [, wavele] [, params] [, parepl] $ ; [, pscale=pscale] [, tausil=tausil] $ ; [, no_uvp=no_uvp] [, fluxes=fluxes] $ ; [, /clumps] [, /reclmp] [, sizclp=sizclp] $ ; [, ranvec=ranvec]) ; ; REQUIRED INPUTS: ; none ; ; OPTIONAL INPUTS: ; arrsiz 2 element array specifying the x and y size of the result ; wavele calculate the model for these wavelengths in micron ; params 16 element array holding the parameters of the model ; parepl 16 element array indicating which parameters are coupled ; ; KEYWORDS: ; pscale pixel scale in milliarcsec (mas) per pixel ; tausil structure holding the silicate optical depth ; no_uvp do not calculate the UV plane, i.e. skip Fourier transform ; clumps make the second component clumpy ; reclmp use provided RANVEC to create the specific clumpiness ; sizclp size of clumps in milliarcseconds (mas) ; ranvec random vector to generate clumpy second component ; ; OUTPUTS: ; gmodel model struct with the 2 component black body model ; fluxes integrated fluxes of the two components ; ranvec random vector used to generate clumpy second component ; ; DESCRIPTION AND EXAMPLE: ; This function calculates a model consisting of two elliptical Gaussian ; black body emitters for several wavelengths. For each of the two ; emitters the following 8 parameters can be specified: offset in RA ; and DEC (in mas), full width at half maximum of the major axis, in- ; clination angle of the ellipsoid. position angle of the ellipsoid ; (counted counterclockwise), optical depth of the silicate absorption, ; temperature and scaling factor. The total number of 16 parameters are ; provided through PARAMS. The wavelengths for which the model is calcu- ; lated are given by the parameter WAVELE. The model can take into ; account silicate obsorption for which the optical depth has to be ; provided through the parameter TAUSIL in form of a corresponding ; structure. Additionally a "screen of clumpiness" can be computed ; simulating a clumpy or inomogeneous flux distribution for the larger ; Gaussian component. The model is returned as a model structure ; compliant to a FITS structure. The values in the image plane are in ; Jy per pixel. Unless the keyword NO_UVP is set, also the UV plane for ; the model is calculated. Using the default settings, the function is ; called by: ; gmodel = midi_2bbmdl() ; ; CALLED BY: ; midi_2bbdat ; midi_2bbcmp ; ; CALLING: ; Makes use of routines for FITS headers from the ASTROLIB library. ; modl_egauss ; modl_mkhead ; modl_m2visi ; ; MODIFICATION HISTORY: ; 2007-02-07 Adapted from earlier test_5.pro by Konrad R. Tristram ; 2009-07-22 K. Tristram: Renamed and adapted to the new fitting routine. ; FUNCTION MIDI_2BBMDL, ARRSIZ, WAVELE, PARAMS, PAREPL, PSCALE=PSCALE, $ TAUSIL=TAUSIL, FLUXES=FLUXES, NO_UVP=NO_UVP, $ CLUMPS=CLUMPS, RECLMP=RECLMP, SIZCLP=SIZCLP, $ RANVEC=RANVEC ; DEFINE A FEW PHYSICAL CONSTANTS ; ------------------------------------------------------------------------------ h_const = 6.62606800d-34 ; Planck's constant c_const = 2.99792458d+08 ; speed of light k_const = 1.38065000d-23 ; Boltzmann's constant ; CHECK IF ARRAY SIZE, WAVELENGTHS AND PIXELSCALE ARE SET ; ------------------------------------------------------------------------------ if (n_elements(arrsiz) eq 0) then arrsiz = [1, 1] * 512 else $ if (n_elements(arrsiz) eq 1) then arrsiz = [1, 1] * arrsiz if (n_elements(wavele) eq 0) then wavele = findgen(21)/4. + 8.0 if (n_elements(pscale) eq 0) then pscale = 2.0 ; Note that the values should have the following contraints: ; pscale < 8 mas / px, otherwise small components are not sampled correctly ; pscale*arrsize > 256 mas, otherwise flux gets lost and UV resolution is bad ; CHECK IF THE CORRECT NUMBER OF PARAMETERS IS GIVEN ; ------------------------------------------------------------------------------ if (n_elements(params) ne 16) then begin print, 'Parameter array must have 16 elements. Using default values.' params = [ 0.0, 0.0, 20.0, 50.0, 40.0, 1.0, 400.0, 1.0, $ 10.0, -5.0, 70.0, 60.0, 130.0, 2.0, 300.0, 0.2] endif ; REPLACE THE PARAMETERS WHICH ARE COUPLED ; ------------------------------------------------------------------------------ if (n_elements(parepl) ne 16) then parepl = indgen(16) parnew = params[parepl] ; START GENERATING A DEFAULT HEADER FOR THE IMAGES ;------------------------------------------------------------------------------- header = modl_mkhead(arrsiz, 'MIDI_2BBMDL', '2 Black Body Gaussians', pscale) ; ADD INPUT PARAMETERS AS HEADER KEYWORDS ;------------------------------------------------------------------------------- str = ' Model parameter component 1: ' sxaddpar, header, 'XOFFS1', parnew [00], str+'offset in RA', before='COMMENT' sxaddpar, header, 'YOFFS1', parnew [01], str+'offset in DEC', before='COMMENT' sxaddpar, header, 'FWHMA1', parnew [02], str+'FWHM', before='COMMENT' sxaddpar, header, 'IANGL1', parnew [03], str+'incl. angle', before='COMMENT' sxaddpar, header, 'PANGL1', parnew [04], str+'position angle', before='COMMENT' sxaddpar, header, 'SILDE1', parnew [05], str+'silicate depth', before='COMMENT' sxaddpar, header, 'TEMPE1', parnew [06], str+'temperature', before='COMMENT' sxaddpar, header, 'FACTO1', parnew [07], str+'scaling factor', before='COMMENT' str = ' Model parameter component 2: ' sxaddpar, header, 'XOFFS2', parnew [08], str+'offset in RA', before='COMMENT' sxaddpar, header, 'YOFFS2', parnew [09], str+'offset in DEC', before='COMMENT' sxaddpar, header, 'FWHMA2', parnew [10], str+'FWHM', before='COMMENT' sxaddpar, header, 'IANGL2', parnew [11], str+'incl. angle', before='COMMENT' sxaddpar, header, 'PANGL2', parnew [12], str+'position angle', before='COMMENT' sxaddpar, header, 'SILDE2', parnew [13], str+'silicate depth', before='COMMENT' sxaddpar, header, 'TEMPE2', parnew [14], str+'temperature', before='COMMENT' sxaddpar, header, 'FACTO2', parnew [15], str+'scaling factor', before='COMMENT' ; CREATE THE STRUCTURE HOLDING THE RESULT ;------------------------------------------------------------------------------- result = {object:'2 Black Body Gaussians', wavele:0.0, arrsiz:arrsiz, $ pscale:pscale, header:header, imgdat:fltarr(arrsiz)} result = replicate(result, n_elements(wavele)) ; COPY THE WAVELENGTHS INTO THE RESULT STRUCTURE ; ------------------------------------------------------------------------------ result.wavele = wavele nu = c_const/(wavele * 1.0e-6) ; CREATE THE ARRAY CONTAINING THE TOTAL FLUXES OF THE INDIVIDUAL COMPONENTS ; ------------------------------------------------------------------------------ fluxes = make_array(2, n_elements(wavele), value=0.0, /float) ; SET THE SIZE OF CLUMPS AND CALCULATE SIZE OF ARRAY FOR CLUMP GENERATION ; ------------------------------------------------------------------------------ if n_elements(sizclp) lt 1 then sizclp = 10. ; size of clumps in mas carrsiz = pscale * arrsiz / sizclp ; CREATE A RANDOM VECTOR IF GENERATING A NEW DISTRIBUTION OF CLUMPS ; ------------------------------------------------------------------------------ if keyword_set(clumps) then begin ranvec = randomu(seed, carrsiz * parnew[15]) reclmp = 1 endif ; CONVERT RANDOM VECTOR INTO A "SCREEN OF CLUMPS" ; ------------------------------------------------------------------------------ if keyword_set(reclmp) then begin clmparr = make_array(carrsiz, carrsiz, /float, value=0) for i=0L, n_elements(ranvec)-1 do begin $ allowed = where(clmparr eq 0) &$ clmparr[allowed[ranvec[i]*n_elements(allowed)]] = 1 &$ endfor kernel = modl_egauss(10, 4., 0, 0, /notime) kernel = kernel/total(kernel) clmparr = convol(clmparr, kernel, /center, /edge_trunc) clmparr = congrid(clmparr, arrsize, arrsize, /interp, cubic=1) clmparr = clmparr > 0 endif ; CALCULATE THE TWO ELLIPTICAL GAUSSIAN DISTRIBUTIONS ; ------------------------------------------------------------------------------ e1 = modl_egauss(arrsiz, parnew[02]/pscale, parnew[03], parnew[04], $ xoffse=parnew[00]/pscale, yoffse=parnew[01]/pscale, /silent) e2 = modl_egauss(arrsiz, parnew[10]/pscale, parnew[11], parnew[12], $ xoffse=parnew[08]/pscale, yoffse=parnew[09]/pscale, /silent) ; DETERMINE THE OPTICAL DEPTH ; ------------------------------------------------------------------------------ if n_elements(tausil) ne 0 then optdepth = interpol(tausil.optdepth, $ tausil.wave, 1e-6 * wavele) > 0.0 else optdepth = 0 ; CALCULATE FLUX OF FIRST COMPONENT ; ------------------------------------------------------------------------------ f1 = 2 * h_const * nu^3 / c_const^2 / $ (exp(h_const * nu / k_const / parnew[06]) - 1) f1 = f1 * (2.0 * !dpi / (3600. * 360. * 1000.))^2 * 1.0e26 f1 = parnew[07] * f1 * exp(- optdepth * parnew[05]) * pscale^2 fluxes[0,*] = f1 * total(e1) ; CALCULATE FLUX OF SECOND COMPONENT ; ------------------------------------------------------------------------------ f2 = 2 * h_const * nu^3 / c_const^2 / $ (exp(h_const * nu / k_const / parnew[14]) - 1) f2 = f2 * (2.0 * !dpi / (3600.*360.*1000.))^2 * 1.0e26 f2 = f2 * exp(- optdepth * parnew[13]) * pscale^2 if keyword_set(reclmp) then e2 = clmparr * e2 else e2 = parnew[15] * e2 fluxes[1,*] = f2 * total(e2) ; LOOP OVER WAVELENGTHS TO MAKE THE FINAL FLUX DISTRIBUTION ; ------------------------------------------------------------------------------ for i=0L,n_elements(wavele)-1 do begin ; COPY THE DATA ARRAY TO THE RESULT STRUCT ;----------------------------------------------------------------------- result[i].imgdat = f1[i] * e1 + f2[i] * e2 ; UPDATE THE WAVELENGTH INFORMATION IN THE HEADER ;----------------------------------------------------------------------- sxaddpar, header, 'WAVELE', wavele[i] sxaddpar, header, 'LAMBDA1', wavele[i] * 1.0E-6 ; COPY THE UPDATED HEADER TO THE RESULT STRUCT ;----------------------------------------------------------------------- result[i].header = header endfor ; RETURN THE MODEL STRUCT WITH OR WITHOUT THE UV PLANE ; ------------------------------------------------------------------------------ if keyword_set(no_uvp) then return, result $ else return, modl_m2visi(result, /silent) END ;+ ; NAME: ; MIDI_2BBDAT ; ; PURPOSE: ; Calculate and return the "visibility data" for the model. ; ; CATEGORY: ; MIDI two component black body fit. ; ; CALLING SEQUENCE: ; mdlvis = midi_2bbdat(transf, params [, /clumps] [, /reclmp] $ ; [, sizclp=sizclp] [, ranvec=ranvec]) ; ; REQUIRED INPUTS: ; transf transfer variable containing information on data and model ; params array holding the parameters of the model ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; clumps make second component clumpy ; reclmp use provided RANVEC to create the specific clumpiness ; sizclp size of clumps in milliarcseconds (mas) ; ranvec random vector to generate clumpy second component ; ; OUTPUTS: ; mdlvis "visibility data" as a simple array ; ; DESCRIPTION AND EXAMPLE: ; This function calculates a two component Gaussian black body model ; using the function MIDI_2BBMDL and returns "visibility data" for it ; in a simple array. The function is mainly intended for being called ; by the fitting routine MPFITFUN. The parameters for the model are ; transferred through the variable PARAMS, all additional information ; needed for the generation of the model and for the extraction of the ; "visibility data" are transferred using the parameter TRANSF. There- ; fore, TRANSF must contain at least (1) the wavelengths for which the ; fit is to be carried out (TRANSF.WAVELE), (2) a three element array ; specifying what data is to be used for the fitting (TRANSF.SELECT) ; and (3) a structure containing the measured data as a MIDI result ; structure (TRANSF.DATVIS). Optionally clumpiness can be simulated. ; Depending on the values in TRANSF.SELECT the "visibility data" re- ; turned contains one or more of (1) the correlated fluxes, (2) the ; visibilities and (3) the differential phases. The function is called ; by: 'fluxes = fit2cbbvisi(transf, params)'. ; ; CALLED BY: ; midi_2bbfit ; ; CALLING: ; midi_2bbmdl ; midi_extvis ; midi_dphase ; ; MODIFICATION HISTORY: ; 2007-02-07 Adapted from earlier test_5.pro by Konrad R. Tristram ; 2009-07-22 K. Tristram: Renamed and adapted to the new fitting routine. ; FUNCTION MIDI_2BBDAT, TRANSF, PARAMS, CLUMPS=CLUMPS, RECLMP=RECLMP, $ SIZCLP=SIZCLP, RANVEC=RANVEC ; IF TRANSF CONTAINS THE CORRECT FIELDS THEN USE THE INFORMATION FOR THE FIT ; ------------------------------------------------------------------------------ if total(tag_names(transf) eq 'ARRSIZ') then arrsiz = transf.arrsiz if total(tag_names(transf) eq 'PSCALE') then pscale = transf.pscale if total(tag_names(transf) eq 'TAUSIL') then tausil = transf.tausil ; CALCULATE THE MODEL ; ------------------------------------------------------------------------------ uvplan = midi_2bbmdl(arrsiz, transf.wavele, params, pscale=pscale, $ tausil=tausil, sizclp=sizclp, clumps=clumps, $ reclmp=reclmp, ranvec=ranvec) ; EXTRACT THE VISIBILITIES AND REMOVE FIRST AND SECOND ORDER PHASES ; ------------------------------------------------------------------------------ mdlvis = midi_dphase(midi_extvis(uvplan, transf.datvis), /silent) ; CONSTRUCT THE RESULT ARRAY ACCORDING TO THE SELECTION IN TRANSF.SELECT ; ------------------------------------------------------------------------------ result = make_array(1, n_elements(mdlvis)) if transf.select[0] then result = [result, mdlvis.corflx] if transf.select[1] then result = [result, mdlvis.visamp] if transf.select[2] then result = [result, mdlvis.visphi] ; RETURN THE RESULT ARRAY ; ------------------------------------------------------------------------------ return, result[1:*, *] END ;+ ; NAME: ; MIDI_2BBFIT ; ; PURPOSE: ; Run the fit for the 2 component Gaussian black body model. ; ; CATEGORY: ; MIDI two component black body fit. ; ; CALLING SEQUENCE: ; midi_2bbfit, transf [, result] ; ; REQUIRED INPUTS: ; transf transfer variable containing information on data and model ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; none ; ; OUTPUTS: ; result array holding the parameters of the best fit model ; ; DESCRIPTION AND EXAMPLE: ; This programme carries out a Levenberg-Marquardt least-squares fit ; of a two component black body emitter to MIDI data. All quantities ; necessary for the fit must be provided using the transfer variable ; TRANSF, which can be created for example by 'transf = midi_2bbpre()' ; (see also there). TRANSF must contain at least (1) the wavelengths for ; which the fit is to be carried out (TRANSF.WAVELE), (2) a three ; element array specifying what data is to be used for the fitting ; (TRANSF.SELECT) and (3) a structure containing the measured data as ; a MIDI result structure (TRANSF.DATVIS). The programme first initia- ; lises the parameter information, then constructs a comparison array ; from the measured data and then carries out the fitting using MPFITFUN ; together with the function MIDI_2BBDAT. The best fit parameters are ; printed out and returned in RESULT. The programme rquires the function ; MPFITFUN to be present. To run the programme call: ; transf = midi_2bbpre() ; midi_2bbfit, transf, result ; ; CALLED BY: ; none ; ; CALLING: ; mpfitfun ; midi_2bbdat ; ; MODIFICATION HISTORY: ; 2007-02-07 Adapted from earlier test_5.pro by Konrad R. Tristram ; 2009-07-22 K. Tristram: Renamed and rewritten completely. ; PRO MIDI_2BBFIT, TRANSF, RESULT ; SET TIMER FOR CALCULATION OF FIT DURATION ; ------------------------------------------------------------------------------ time = systime(/seconds) ; CHECK IF TRANSF IS A STRUCTURE ; ------------------------------------------------------------------------------ if size(transf, /type) ne 8 then begin print, 'TRANSF must be a structure containing specifications for ' + $ 'the fit. Returning.' return endif ; INITIALISE THE PARAMETERS TO BE FITTED ; ------------------------------------------------------------------------------ parinfo = replicate({value:0.d0, fixed:0, limited:[1,1], relstep:0.1d0, $ mpside:2, step:0, limits:[0.D,0]}, 16) parinfo[00].value = 0.0 ;x-offset1 parinfo[00].limits[0] = -100.0 parinfo[00].limits[1] = 100.0 parinfo[01].value = 0.0 ;y-offset1 parinfo[01].limits[0] = -100.0 parinfo[01].limits[1] = 100.0 parinfo[02].value = 20.0 ;fwhm1 parinfo[02].limits[0] = 0.0 parinfo[02].limits[1] = 500.0 parinfo[03].value = 60.0 ;iangle1 parinfo[03].limits[0] = 0.0 parinfo[03].limits[1] = 90.0 parinfo[04].value = 60.0 ;pangle1 parinfo[04].limits[0] = -270.0 parinfo[04].limits[1] = 270.0 parinfo[05].value = 1.0 ;sildep1 parinfo[05].limits[0] = -10.0 parinfo[05].limits[1] = 10.0 parinfo[06].value = 400.0 ;temp1 parinfo[06].limits[0] = 0.0 parinfo[06].limits[1] = 1500.0 parinfo[07].value = 1.0 ;fact1 parinfo[07].limits[0] = 0.0 parinfo[07].limits[1] = 1.0 parinfo[08].value = 0.0 ;x-offset2 parinfo[08].limits[0] = -100.0 parinfo[08].limits[1] = 100.0 parinfo[09].value = 0.0 ;y-offset2 parinfo[09].limits[0] = -100.0 parinfo[09].limits[1] = 100.0 parinfo[10].value = 90.0 ;fwhm2 parinfo[10].limits[0] = 0.0 parinfo[10].limits[1] = 500.0 parinfo[11].value = 15.0 ;iangle2 parinfo[11].limits[0] = 0.0 parinfo[11].limits[1] = 90.0 parinfo[12].value = 60.0 ;pangle2 parinfo[12].limits[0] = -270.0 parinfo[12].limits[1] = 270.0 parinfo[13].value = 1.0 ;sildep2 parinfo[13].limits[0] = -10.0 parinfo[13].limits[1] = 10.0 parinfo[14].value = 300.0 ;temp2 parinfo[14].limits[0] = 0.0 parinfo[14].limits[1] = 1500.0 parinfo[15].value = 1.0 ;fact2 parinfo[15].limits[0] = 0.0 parinfo[15].limits[1] = 1.0 if total(tag_names(transf) eq 'PARAMS') then $ if n_elements(transf.params) eq 16 then parinfo.value = transf.params if total(tag_names(transf) eq 'PARFIX') then $ if n_elements(transf.parfix) eq 16 then parinfo.fixed = transf.parfix ; CONSTRUCT THE MEASUREMENT AND ERROR ARRAYS ACCORDING TO THE SELECTION ; ------------------------------------------------------------------------------ datcom = midi_datcom(transf.datvis, transf.wavele, 5) datvis = make_array(1, n_elements(transf.datvis)) daterr = datvis if transf.select[0] then begin datvis = [datvis, datcom.corflx] daterr = [daterr, datcom.corflxerr] endif if transf.select[1] then begin datvis = [datvis, datcom.visamp] daterr = [daterr, datcom.visamperr] endif if transf.select[2] then begin datvis = [datvis, datcom.visphi] daterr = [daterr, datcom.visphierr] endif datvis = datvis[1:*, *] daterr = daterr[1:*, *] ; RUN THE FIT USING MPFITFUN ; ------------------------------------------------------------------------------ result = mpfitfun('midi_2bbdat', transf, datvis, daterr, parinfo.value, $ parinfo=parinfo, xtol=1.d-20, status=status, $ bestnorm=bestnorm, perror=perror, maxiter=20) ; PRINT OUT THE RESULT ; ------------------------------------------------------------------------------ dof = (n_elements(datvis) - n_elements(where(parinfo.fixed eq 0))) print, ' value fitted' print, ' ----- ------' varnam = ['x-offset1', 'y-offset1', 'fwhm1', 'iangle1', 'pangle1', $ 'sildep1', 'temp1', 'fact1', $ 'x-offset2', 'y-offset2', 'fwhm2', 'iangle2', 'pangle2', $ 'sildep2', 'temp2', 'fact2'] for i=0,n_elements(result)-1 do begin print, varnam[i]+':', result[i], perror[i], parinfo[i].fixed, $ format='(A-13, F10.3, " ±", F10.3, I6)' endfor print, bestnorm , format='("chi-squared: ", F10.3)' print, bestnorm / dof , format='("reduced X^2: ", F10.3)' print, dof , format='("dof.: ", I6)' print, 'Computation time: ', systime(/SECONDS) - time, ' sec' ;print, result, format='(16(G7.5, :, ","))' END ;+ ; NAME: ; MIDI_2BBCMP ; ; PURPOSE: ; Compare a 2 component Gaussian black body model to the data. ; ; CATEGORY: ; MIDI two component black body fit. ; ; CALLING SEQUENCE: ; midi_2bbcmp, transf, params [, wavele] [, wavran=wavran] [, /corflx] ; [, /visamp][, /visphi] [, /normal] [, pmulti=pmulti] ; [, psfile=psfile] [, /compact] ; ; REQUIRED INPUTS: ; transf transfer variable containing information on data and model ; params array holding the parameters of the model ; ; OPTIONAL INPUTS: ; wavele calculate the model for these wavelengths in micron ; ; KEYWORDS: ; wavran the wavelength range for which the plots will be produced ; corflx make comparison plots for the correlated fluxes ; visamp make comparison plots for the visibilities ; visphi make comparison plots for the differential phases ; normal normalise the scaling of the y axis in all plots ; pmulti number of colunms and rows per page in the plot ; psfile prefix of the PS file for POSTSCRIPT output ; compact make a compact plot in the case of PS output ; ; OUTPUTS: ; - plots comparing the measured data to the model data ; ; DESCRIPTION AND EXAMPLE: ; This programme calculates a 2 component Gaussian black body model ; according to the parameters in PARAMS and the settings specified by ; WAVELE and/or TRANSF. Two additional models are calculated for each ; of the two components alone. Then the programme MIDI_PLTCMP is called ; to carry out the comparison of the model to the data. The programme ; is simply called by: ; midi_2bbcmp, transf, params ; ; CALLED BY: ; none ; ; CALLING: ; midi_2bbmdl ; midi_extvis ; midi_dphase ; midi_pltcmp ; ; MODIFICATION HISTORY: ; 2009-08-03 Written by Konrad R. W. Tristram ; 2009-08-18 K. Tristram: Added keyword PMULTI. ; PRO MIDI_2BBCMP, TRANSF, PARAMS, WAVELE, WAVRAN=WAVRAN, CORFLX=CORFLX, $ VISAMP=VISAMP, VISPHI=VISPHI, NORMAL=NORMAL, PMULTI=PMULTI, $ PSFILE=PSFILE, COMPACT=COMPACT ; SELECT THE WAVELENGTH TO CALCULATE THE MODEL FOR ; ------------------------------------------------------------------------------ if n_elements(wavele) lt 2 then wavele=transf.wavele ; CALCULATE THE FULL MODEL ; ------------------------------------------------------------------------------ uvplan = midi_2bbmdl(transf.arrsiz, wavele, params, $ pscale=transf.pscale, tausil=transf.tausil) visisn = midi_dphase(midi_extvis(uvplan, transf.datvis)) ; CALCULATE THE MODEL FOR THE SECOND COMPONENT ONLY ; ------------------------------------------------------------------------------ param1 = params param1[[15]] = 0.0 uvpla1 = midi_2bbmdl(transf.arrsiz, wavele, param1, $ pscale=transf.pscale, tausil=transf.tausil) visis1 = midi_dphase(midi_extvis(uvpla1, transf.datvis)) ; CALCULATE THE MODEL FOR THE FIRST COMPONENT ONLY ; ------------------------------------------------------------------------------ param2 = params param2[[07]] = 0.0 uvpla2 = midi_2bbmdl(transf.arrsiz, wavele, param2, $ pscale=transf.pscale, tausil=transf.tausil) visis2 = midi_dphase(midi_extvis(uvpla2, transf.datvis)) ; NO DO THE ACTUAL COMPARISON FORWARDING ALL KEYWORDS FOR MIDI_PLTCMP ; ------------------------------------------------------------------------------ midi_pltcmp, transf.datvis, visisn, visis1, visis2, wavran=wavran, $ corflx=corflx, visamp=visamp, visphi=visphi, normal=normal, $ pmulti=pmulti, psfile=psfile, compact=compact END