; MODL_EGAUSS - Generate a generalised elliptical Gaussian or uniform disk. ; result = modl_egauss(arrsiz, fwhmax, iangle, pangle, /ufdisk) ; ; MODL_GAUSSM - Generate an set of generalised elliptical Gaussians. ; result = modl_gaussm(arrsiz, wavele, fwhmax, iangle, pangle) ; ; MODL_POLARM - Generate a polar model for a single wavelength. ; result = modl_polarm(arrsiz, fwhmax, iangle, pangle, /gaussd) ; ; MODL_TDISKM - Generate a temperature disk model with T(r) = r^(-EXPONE) ; result = modl_tdiskm(arrsiz, wavele, innrad, iangle, pangle) ; ;+ ; NAME: ; MODL_EGAUSS ; ; PURPOSE: ; Generate a generalised elliptical Gaussian or uniform disk. ; ; CATEGORY: ; MIDI models. ; ; CALLING SEQUENCE: ; result = modl_egauss([arrsiz] [, fwhmax] [, iangle] [, pangle] ; [, boxyns=boxyns] [, xoffse=xoffse] ; [, yoffse=yoffse] [, expone=expone] ; [, /ufdisk] [, /silent]) ; ; REQUIRED INPUTS: ; None. ; ; OPTIONAL INPUTS: ; arrsiz 2 element array specifying the x and y size of the result ; fwhmax FWHM of the longer axis of the elliptical Gaussian ; iangle inclination angle of the model, governs ellipticity ; pangle position angle of the ellipsoid counting counter clockwise ; ; KEYWORDS: ; boxyns boxyness parameter for boxy or disky shape ; xoffse offset from the centre of the array in x direction ; yoffse offset from the centre of the array in y direction ; expone exponent for the radial dependency in the EXP function ; ufdisk calculate a uniform disk instead of a Gaussian ; silent suppress the output of the time the programme took to run ; ; OUTPUTS: ; result the resulting model as two dimensional array ; ; DESCRIPTION AND EXAMPLE: ; This function creates a generalised elliptical Gaussian brightness ; distribution including boxy and disky profiles as well as uniform ; disks. The size, ellipticity and orientation of the Gaussian distri- ; bution are governed by the parameters FWHMAX, IANGLE and PANGLE, ; respectively. For BOXYNS < 0 disky distributions, for BOXYNS = 0 a ; normal elliptical distribution and for BOXYNS > 0 boxy distributions ; are obtained. For BOXYNS > 100 a truely boxy distribution is calcu- ; lated. The position of the brightness distribution can be offset from ; the centre of the array by specifying XOFFSE and YOFFSE. By varying ; EXPONE more bell shaped or more centrally concentrated distributions ; can be obtained. For EXPONE -> Inf the result becomes a uniform disk, ; which can directly be obtained by using the keyword UFDISK. The ; simplest way to call the function is by invocing 'result = ellgauss()'. ; ; CALLED BY: ; modl_gaussm ; midi_2bbmdl ; ; CALLING: ; none ; ; MODIFICATION HISTORY: ; 2005-09-04 Written by Konrad R. W. Tristram ; 2006-01-12 K. Tristram: Removed error in FWHM to SIGMA conversion. ; 2006-04-28 K. Tristram: Added support of boxyness and diskyness. ; 2006-06-19 K. Tristram: Added possibility for flat disk distribution. ; 2009-07-13 K. Tristram: Renamed function and various other changes. ; FUNCTION MODL_EGAUSS, ARRSIZ, FWHMAX, IANGLE, PANGLE, BOXYNS=BOXYNS, $ XOFFSE=XOFFSE, YOFFSE=YOFFSE, EXPONE=EXPONE, $ UFDISK=UFDISK, SILENT=SILENT ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- systim = systime(/seconds) 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(fwhmax) eq 0) then fwhmax = 20.0 if (n_elements(iangle) eq 0) then iangle = 60.0 if (n_elements(pangle) eq 0) then pangle = 30.0 if (n_elements(boxyns) eq 0) then boxyns = 0.0 if (n_elements(xoffse) eq 0) then xoffse = 0.0 if (n_elements(yoffse) eq 0) then yoffse = 0.0 if (n_elements(expone) eq 0) then expone = 2.0 boxyn2 = float(boxyns) + 2.0 ; CALCULATE WIDTH OF MINOR AND MAJOR AXIS OF ELLIPSOID ;------------------------------------------------------------------------------- if keyword_set(ufdisk) then gsigma = float(fwhmax) $ else gsigma = float(fwhmax) / sqrt(8.0 * alog(2.0)) x_sigm = gsigma * cos(iangle/180.0*!pi) y_sigm = gsigma ; GENERATE RADIAL DEPENDENCY TERM AND MAKE IT ELLIPTICAL ;------------------------------------------------------------------------------- x = findgen(arrsiz[0]) y = findgen(arrsiz[1]) x = (x-arrsiz[0]/2-xoffse) # replicate(1.0, arrsiz[1]) y = replicate(1.0, arrsiz[0]) # (y-arrsiz[1]/2-yoffse) ; ADD THE POSITION ANGLE DEPENDENCY ;------------------------------------------------------------------------------- s = sin(pangle/180.0*!pi) & c = cos(pangle/180.0*!pi) t = ( x * c + y * s) / x_sigm y = (-x * s + y * c) / y_sigm x = temporary(t) ; CALCULATE RESULTING ARRAY, I.E. EXPONENTIAL PART OF GAUSSIAN OR DISK AREA ;------------------------------------------------------------------------------- if keyword_set(ufdisk) then begin result = make_array(arrsiz, value=0.0*fwhmax) index = where((abs(x)^boxyn2 + abs(y)^boxyn2)^(1./boxyn2) le 0.5) if index[0] ge 0 then result[index] = 1.0 endif else begin ; FOR VERY LARGE BOXYNS ASSUME THAT REAL BOX SHAPE IS DESIRED ;----------------------------------------------------------------------- if boxyn2 ge 102 then begin result = exp(-0.5 * (abs(x) > abs(y))^expone) ; OTHERWISE CALCULATE GENERALISED ELLIPTICAL GAUSSIAN ;----------------------------------------------------------------------- endif else begin result = exp(-0.5*(abs(x)^boxyn2+abs(y)^boxyn2)^(expone/boxyn2)) endelse endelse ; PRINT TIME PROGRAMME TOOK TO RUN ;------------------------------------------------------------------------------- if not keyword_set(silent) then begin print, 'Computation time: ', systime(/seconds) - systim, ' sec' endif ; RETURN RESULT ;------------------------------------------------------------------------------- return, result END ;+ ; NAME: ; MODL_GAUSSM ; ; PURPOSE: ; Generate an set of generalised elliptical Gaussians. ; ; CATEGORY: ; MIDI models. ; ; CALLING SEQUENCE: ; result = modl_gaussm([, arrsiz] [, wavele] [, fwhmax] [, iangle] $ ; [, pangle] [, boxyns=boxyns] [, pscale=pscale] $ ; [, xoffse=xoffse] [, yoffse=yoffse] [, /ufdisk] $ ; [, expone=expone] [, /silent] [, /output]) ; ; 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 ; fwhmax FWHM of the longer axis of the elliptical Gaussian ; iangle inclination angle of the model, governs ellipticity ; pangle position angle of the ellipsoid counting counter clockwise ; ; KEYWORDS: ; boxyns boxyness parameter for boxy or disky shape ; pscale pixel scale: milliarcsec (mas) per pixel ; xoffse offset from the centre of the array in x direction ; yoffse offset from the centre of the array in y direction ; expone exponent for the radial dependency in the EXP function ; ufdisk calculate a uniform disk instead of a Gaussian ; silent suppress the output of the time the programme took to run ; output make plot of the flux distribution for every wavelength ; ; OUTPUTS: ; result resulting model as FITS struct ; ; DESCRIPTION AND EXAMPLE: ; This function creates a generalised elliptical Gaussian model with ; different sizes for different wavelengths. Currently the dependence ; of the size on wavelength is hardcoded: it is a simple linear ; relationship. The function uses the function MODL_EGAUSS to calc- ; ulate the model and all keywords of this function are also accepted ; by MODL_GAUSSM. The simplest way to call the function is by ; 'result = gaussmodel()'. ; ; CALLED BY: ; none ; ; CALLING: ; Makes use of routines for FITS headers from the ASTROLIB library. ; modl_egauss ; modl_mkhead ; modl_m2plot ; ; MODIFICATION HISTORY: ; 2005-09-04 Written by Konrad R. W. Tristram ; 2006-01-12 K. Tristram: Added support for pixelscale, offset and output. ; 2006-05-31 K. Tristram: Added support of boxyness and diskyness. ; 2006-06-26 K. Tristram: Made RA axis count in correct direction in plot. ; 2009-07-13 K. Tristram: Renamed function and various other changes. ; FUNCTION MODL_GAUSSM, ARRSIZ, WAVELE, FWHMAX, IANGLE, PANGLE, BOXYNS=BOXYNS, $ PSCALE=PSCALE, XOFFSE=XOFFSE, YOFFSE=YOFFSE, $ EXPONE=EXPONE, UFDISK=UFDISK, SILENT=SILENT, $ OUTPUT=OUTPUT ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- systim = systime(/seconds) 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 = [8.0, 9.0, 10.0, 11.0, 12.0, 13.0] if (n_elements(fwhmax) eq 0) then fwhmax = 20.0 if (n_elements(pscale) eq 0) then pscale = 1.0 if (n_elements(xoffse) eq 0) then xoffse = 0.0 if (n_elements(yoffse) eq 0) then yoffse = 0.0 if ~ keyword_set(ufdisk) then ufdisk=0B ; START GENERATING A DEFAULT HEADER FOR THE IMAGES ;------------------------------------------------------------------------------- header = modl_mkhead(arrsiz, 'MODL_GAUSSM', 'Gaussian Model', pscale) ; ADD INPUT PARAMETERS AS HEADER KEYWORDS ;------------------------------------------------------------------------------- str = ' Model parameter: ' sxaddpar, header, 'FWHMAX', 0.0, str+'FWHM', before='COMMENT' sxaddpar, header, 'IANGLE', 0.0, str+'inclination angle', before='COMMENT' sxaddpar, header, 'PANGLE', 0.0, str+'position angle', before='COMMENT' sxaddpar, header, 'BOXYNS', 0.0, str+'boxyness', before='COMMENT' sxaddpar, header, 'XOFFSE', xoffse, str+'x offset from centre', before='COMMENT' sxaddpar, header, 'YOFFSE', yoffse, str+'y offset from centre', before='COMMENT' sxaddpar, header, 'EXPONE', 0.0, str+'exponent', before='COMMENT' if keyword_set(ufdisk) then sxaddpar, header, 'UFDISK', 'T', ' Model is a ' + $ 'uniform disk', before='COMMENT' $ else sxaddpar, header, 'UFDISK', 'F', ' Model is ' + $ 'not a uniform disk', before='COMMENT' ; CREATE THE STRUCT HOLDING THE RESULT ;------------------------------------------------------------------------------- result = {object:'Gaussian Model', wavele:0.0, arrsiz:arrsiz, pscale:pscale, $ header:header, imgdat:make_array(arrsiz, value=0.0)} result = replicate(result, n_elements(wavele)) ; LOOP OVER THE WAVELENGTHS ;------------------------------------------------------------------------------- for i=0, n_elements(wavele)-1 do begin ; COPY THE WAVELENGTHS AND CALCULATE THE IMAGE USING MODL_EGAUSS ;----------------------------------------------------------------------- result[i].wavele = wavele[i] ; factor = (3.8 * wavele[i] + 2.0) / 40. / pscale ; factor = 1. / pscale factor = 1. result[i].imgdat = modl_egauss(arrsiz, (factor * fwhmax), iangle, $ pangle, boxyns=boxyns, expone=expone, $ xoffse=xoffse/pscale, ufdisk=ufdisk, $ yoffse=yoffse/pscale, /silent) ; UPDATE THE HEADER KEYWORDS TO THE VALUES USED ;----------------------------------------------------------------------- sxaddpar, header, 'WAVELE', wavele[i] sxaddpar, header, 'LAMBDA1', wavele[i] * 1.0E-6 sxaddpar, header, 'FWHMAX', (factor * fwhmax) sxaddpar, header, 'IANGLE', iangle sxaddpar, header, 'PANGLE', pangle sxaddpar, header, 'BOXYNS', boxyns sxaddpar, header, 'EXPONE', expone ; COPY THE UPDATED HEADER TO THE RESULT STRUCT ;----------------------------------------------------------------------- result[i].header = header endfor ; PRINT TIME PROGRAMME TOOK TO RUN ;------------------------------------------------------------------------------- if not 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_POLARM ; ; PURPOSE: ; Generate a polar model for a single wavelength. ; ; CATEGORY: ; MIDI models. ; ; CALLING SEQUENCE: ; result = modl_polarm([, arrsiz] [, fwhmax] [, iangle] [, pangle] $ ; [, expone=expone] [, porder=porder] $ ; [, xoffse=xoffse] [, yoffse=yoffse] [, /gaussd] $ ; [, pscale=pscale] [, /silent] [, /output]) ; ; REQUIRED INPUTS: ; None. ; ; OPTIONAL INPUTS: ; arrsiz 2 element array specifying the x and y size of the result ; fwhmax FWHM of the longer axis of the intensity distribution ; iangle strength of asymmetry (ellipticity) in radial direction ; pangle position angle of the model counting counter clockwise ; ; KEYWORDS: ; expone exponent for sigma of the radial dependency ; porder order of the angular dependency ; xoffse offset from the centre of the array in x direction ; yoffse offset from the centre of the array in y direction ; gaussd make the distribution Gaussian in radial direction ; pscale pixel scale: milliarcsec (mas) per pixel ; silent suppress the output of the time the programme took to run ; output make plot of the flux distribution for every wavelength ; ; OUTPUTS: ; result resulting model as FITS struct ; ; DESCRIPTION AND EXAMPLE: ; This function creates a top hat or Gaussian (when /GUASSD is set) ; model from a polar radial dependency. A circle or ellipse is ob- ; tained for EXPONE = 2 and PORDER = 1. EXPONE > 2 gives a boxy ; ellipsoid, while EXPONE > 2 leads to a disky ellipsoid. EXPONE > 3 ; leads to star formed distributions. A higher porder leads to a ; repeating structure in the angular direction. The simplest way to ; call this function is by 'result = gaussmodel()'. ; ; CALLED BY: ; none ; ; CALLING: ; Makes use of routines for FITS headers from the ASTROLIB library. ; modl_mkhead ; modl_m2plot ; ; MODIFICATION HISTORY: ; 2006-05-01 Written by Konrad R. W. Tristram ; 2009-07-13 K. Tristram: Renamed function and various other changes. ; FUNCTION MODL_POLARM, ARRSIZ, FWHMAX, IANGLE, PANGLE, EXPONE=EXPONE, $ PORDER=PORDER, XOFFSE=XOFFSE, YOFFSE=YOFFSE, $ GAUSSD=GAUSSD, PSCALE=PSCALE, SILENT=SILENT, $ OUTPUT=OUTPUT ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- systim = systime(/seconds) 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(fwhmax) eq 0) then fwhmax = 20.0 if (n_elements(iangle) eq 0) then iangle = 60.0 if (n_elements(pangle) eq 0) then pangle = 30.0 if (n_elements(expone) eq 0) then expone = 2.0 if (n_elements(porder) eq 0) then porder = 1.0 if (n_elements(xoffse) eq 0) then xoffse = 0.0 if (n_elements(yoffse) eq 0) then yoffse = 0.0 if (n_elements(pscale) eq 0) then pscale = 1.0 ; SET THE DUMMY WAVELENGTH ;------------------------------------------------------------------------------- wavele = 10.0 ; GENERATE THE X YND Y ARRAYS ;------------------------------------------------------------------------------- x = findgen(arrsiz[0]) y = findgen(arrsiz[1]) x = (x-arrsiz[0]/2-xoffse/pscale) # replicate(1.0, arrsiz[1]) y = replicate(1.0, arrsiz[0]) # (y-arrsiz[1]/2-yoffse/pscale) ; CONVERT CARTESIAN COORDINATES TO POLAR COORDINATES ;------------------------------------------------------------------------------- radius = sqrt(x^2. + y^2.) polang = pangle / 180. * !pi - atan(x, y) polang = porder * polang ; CALCULATE SIGMA FOR GAUSSIAN FROM FWHM AND POLAR DEPENDENCY ;------------------------------------------------------------------------------- sigmsq = fwhmax^2 / pscale^2 sigmsq = sigmsq / ((abs(cos(polang)))^expone + $ (abs(sin(polang)))^expone / cos(iangle/180.0*!pi)) ; CALCULATE WIDTH OF TOP HAT OR EXPONENTIAL PART OF GAUSSIAN ;------------------------------------------------------------------------------- if keyword_set(gaussd) then begin imgdat = exp(-0.5 * radius^2 / sigmsq * (8. * alog(2))) endif else begin imgdat = make_array(arrsiz, /float, value=0.0*fwhmax) index = where((radius^2 / sigmsq * 4) le 1, count) if count ne 0 then imgdat(index) = 1.0 endelse ; START GENERATING A DEFAULT HEADER FOR THE IMAGES ;------------------------------------------------------------------------------- header = modl_mkhead(arrsiz, 'MODL_POLARM', 'Polar Model', pscale) ; ADD INPUT PARAMETERS AS HEADER KEYWORDS ;------------------------------------------------------------------------------- str = ' Model parameter: ' sxaddpar, header, 'FWHMAX', fwhmax, str+'FWHM', before='COMMENT' sxaddpar, header, 'IANGLE', iangle, str+'inclination angle', before='COMMENT' sxaddpar, header, 'PANGLE', pangle, str+'position angle', before='COMMENT' sxaddpar, header, 'EXPONE', expone, str+'exponent', before='COMMENT' sxaddpar, header, 'PORDER', porder, str+'angular order', before='COMMENT' sxaddpar, header, 'XOFFSE', xoffse, str+'x offset from centre', before='COMMENT' sxaddpar, header, 'YOFFSE', yoffse, str+'y offset from centre', before='COMMENT' if keyword_set(gaussd) then sxaddpar, header, 'GAUSSD', 'T', ' Model has ' + $ 'Gaussian radial dependency', before='COMMENT' $ else sxaddpar, header, 'GAUSSD', 'F', ' Model has ' + $ 'top hat radial dependency', before='COMMENT' ; CONVERT RESULT TO A STRUCT FOR FURTHER USAGE ;------------------------------------------------------------------------------- result = {object:'Polar Model', wavele:wavele, arrsiz:arrsiz, pscale:pscale, $ header:header, imgdat:imgdat} ; PRINT TIME PROGRAMME TOOK TO RUN ;------------------------------------------------------------------------------- if not 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_TDISKM ; ; PURPOSE: ; Generate a temperature disk model with T(r) = r^(-EXPONE) ; ; CATEGORY: ; MIDI models. ; ; CALLING SEQUENCE: ; result = modl_tdiskm([arrsiz] [, wavele] [, innrad] [, iangle] $ ; [, pangle] [, xoffse=xoffse] [, yoffse=yoffse] $ ; [, expone=expone] [, t_maxi=t_maxi] $ ; [, silabs=silabs] [, sildep=sildep] $ ; [, pscale=pscale] [, /fillit] [, /behave] $ ; [, /output] [, /silent]) ; REQUIRED INPUTS: ; None. ; ; OPTIONAL INPUTS: ; arrsiz 2 element array specifying the x and y size of the result ; wavele calculate model for these wavelengths in micron ; innrad inner radius of disk, where t_maxi is reached ; iangle inclination angle of the disk model ; pangle position angle of the disk counting counter clockwise ; ; KEYWORDS: ; xoffse offset from the centre of the array in x direction ; yoffse offset from the centre of the array in y direction ; expone exponent for the dependency of temperature on radius ; t_maxi maximum temperature at 'sublimation radius' INNRAD ; silabs data struct containing silicate absorption curve ; sildep column density of absorbing dust ; pscale pixel scale: milliarcsec (mas) per pixel ; fillit fill the central hole by emission of t_maxi ; behave use modified temperature decrease that is well behaved ; at the centre and no central hole ; output make plot of the flux distribution for every wavelength ; silent suppress the output of the time the programme took to run ; ; OUTPUTS: ; result resulting model as FITS struct ; ; DESCRIPTION AND EXAMPLE: ; This programme generates the flux intensity for a simple black body ; disk at certain wavelengths (specified by WAVELE). The disk has an ; inner radius INNRAD where the 'dust sublimation temperature' T_MAXI ; is assumed to be reached. From there the temperature drops according ; to an exponential law with the exponent EXPONE: ; T = r ^(-EXPONE). ; If the keyword BEHAVE is set, then a modified function, which is ; "well behaved" at the nucleus (i.e. no central hole), is used: ; T = r ^(-EXPONE) * ((r^4+1)/r^4)(-EXPONE/4) ; The disk inclination angle is governed by IANGLE, the orientation by ; the position angle PANGLE. In the simplest way a model can be genera- ; ted using the default settings by invoking: 'model = tempdisk()'. ; The silicate absorption has to be provided to the functen in the ; form of a structure containing the wavelength and the optical depth ; as follows: silabs = {wave:float, optdepth:float}. ; For Circinus the settings are: ; distance: 4.2 Mpc -> 1'' ~ 20 pc ; rinner = 0.2 pc ~ 10 mas (Bianchi 2001) ; router = 20.0 pc ~ 1000 mas (Bianchi 2001) ; iangle = 65.0 ; pangle = 42.0 (Greenhill 2004) ; result = modl_tdiskm(512, wavele, 10.0, 65.0, 42.0, expone=0.80, $ ; t_maxi= 500.0, silabs=silabs, sildep=0.5) ; ; CALLED BY: ; none ; ; CALLING: ; Makes use of routines for FITS headers from the ASTROLIB library. ; modl_mkhead ; modl_m2plot ; ; MODIFICATION HISTORY: ; 2005-09-04 Written by Konrad R. W. Tristram ; 2005-09-13 K. Tristram: Changed handling of hole and outside. ; 2005-09-13 K. Tristram: Changed position in code, where hole and outside ; are calculated. Add keywords FILLIT and BEHAVE. ; 2006-01-12 K. Tristram: Added support for pixelscale, offset and output. ; 2006-06-26 K. Tristram: Made RA axis count in correct direction in plot. ; 2006-06-26 K. Tristram: Changed behaviour for keyword BEHAVE: replaced ; Gaussian by different function at nucleus. ; 2009-07-13 K. Tristram: Renamed function and various other changes. ; FUNCTION MODL_TDISKM, ARRSIZ, WAVELE, INNRAD, IANGLE, PANGLE, XOFFSE=XOFFSE, $ YOFFSE=YOFFSE, EXPONE=EXPONE, T_MAXI=T_MAXI, $ SILABS=SILABS, SILDEP=SILDEP, PSCALE=PSCALE, $ FILLIT=FILLIT, BEHAVE=BEHAVE, OUTPUT=OUTPUT, $ SILENT=SILENT ; INITIALISATION AND VARIABLE CHECK ;------------------------------------------------------------------------------- systim = systime(/seconds) 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 = 12.0 if (n_elements(innrad) eq 0) then innrad = 10.0 else innrad *= 1.0 if (n_elements(iangle) eq 0) then iangle = 60.0 if (n_elements(pangle) eq 0) then pangle = 30.0 if (n_elements(xoffse) eq 0) then xoffse = 0.0 if (n_elements(yoffse) eq 0) then yoffse = 0.0 if (n_elements(expone) eq 0) then expone = 0.90 if (n_elements(t_maxi) eq 0) then t_maxi = 500.0 else t_maxi *= 1.0 if (n_elements(sildep) eq 0) then sildep = 0.5 else sildep *= 1.0 if (n_elements(fillit) eq 0) then fillit = 0B if (n_elements(behave) eq 0) then behave = 0B if (n_elements(pscale) eq 0) then pscale = 1.0 h_const = 6.626068d-34 ; Planck's constant c_const = 2.99792458d8 ; speed of light k_const = 1.380650d-23 ; Boltzmann's constant t_mini = 3.0 ; minimum temper ; CALCULATE WIDTH OF MINOR AND MAJOR AXIS OF THE ELLIPSOID ;------------------------------------------------------------------------------- x_sigm = cos(iangle/180.0*!pi) y_sigm = 1.0 ; GENERATE RADIAL DEPENDENCY TERM AND MAKE IT ELLIPTICAL ;------------------------------------------------------------------------------- x = findgen(arrsiz[0]) y = findgen(arrsiz[1]) x = (x-arrsiz[0]/2-xoffse) # replicate(1.0, arrsiz[1]) y = replicate(1.0, arrsiz[0]) # (y-arrsiz[1]/2-yoffse) ; ADD THE POSITION ANGLE DEPENDENCY ;------------------------------------------------------------------------------- s = sin(pangle/180.0*!pi) & c = cos(pangle/180.0*!pi) t = ( x * c + y * s) / x_sigm y = (-x * s + y * c) / y_sigm x = temporary(t) ; CALCULATE TEMPERATURE OF THE DISK AT GIVEN POINT AND DEFINE CENTRAL HOLE ;------------------------------------------------------------------------------- radius = sqrt(x^2 + y^2) / innrad singul = where(radius eq 0.0) if behave then begin ; temper = t_maxi * exp(innrad^expone - (sqrt(x^2 + y^2))^expone) ; !! THIS WAS THE OLD FUNCTION USED FOR THE BEHAVE CASE !! temper = t_maxi * radius^(-expone) * $ ((radius^4 + 1) / radius^4)^(-expone/4) if (singul[0] ge 0) then temper[singul] = t_maxi cntrhole = [-1] endif else begin temper = t_maxi * radius^(-expone) if (singul[0] ge 0) then temper[singul] = !values.f_infinity cntrhole = where(temper gt t_maxi) endelse ; REPLACE VALUES ON THE OUTSIDE BY MINIMUM TEMPERATURE ;------------------------------------------------------------------------------- outsid = where(temper lt t_mini) if (outsid[0] ge 0) then begin temper[outsid] = replicate(t_mini, n_elements(outsid)) endif ; REPLACE VALUES IN CENTRAL HOLE EITHER BY MINIMUM OR BY MAXIMUM TEMPERATURE ;------------------------------------------------------------------------------- if (cntrhole[0] gt 0) then begin if fillit then temper[cntrhole] = t_maxi else temper[cntrhole] = t_mini endif ; TRANSFORM WAVELENGTH INTO FREQUENCY ;------------------------------------------------------------------------------- nu = c_const/(wavele * 1.0e-6) ; START GENERATING A DEFAULT HEADER FOR THE IMAGES ;------------------------------------------------------------------------------- header = modl_mkhead(arrsiz, 'MODL_TDISKM', 'Temperature Disk Model', pscale) ; ADD INPUT PARAMETERS AS HEADER KEYWORDS ;------------------------------------------------------------------------------- str = ' Model parameter: ' sxaddpar, header, 'INNRAD', innrad, str+'inner radius', before='COMMENT' sxaddpar, header, 'IANGLE', iangle, str+'inclination angle', before='COMMENT' sxaddpar, header, 'PANGLE', pangle, str+'position angle', before='COMMENT' sxaddpar, header, 'XOFFSE', xoffse, str+'x offset from centre', before='COMMENT' sxaddpar, header, 'YOFFSE', yoffse, str+'y offset from centre', before='COMMENT' sxaddpar, header, 'EXPONE', expone, str+'exponent', before='COMMENT' sxaddpar, header, 'T_MAXI', t_maxi, str+'temperature at INNRAD', befor='COMMENT' sxaddpar, header, 'SILDEP', sildep, str+'optical depth', before='COMMENT' ; CREATE THE STRUCT HOLDING THE RESULT ;------------------------------------------------------------------------------- result = {object:'Temperature Disk Model', wavele:0.0, arrsiz:arrsiz, $ pscale:pscale, header:header, imgdat:make_array(arrsiz, value=0.0)} result = replicate(result, n_elements(nu)) ; IF SILICATE ABSORPTION IF GIVEN THEN CALCULATE THE ABSORPTION FACTOR ;------------------------------------------------------------------------------- if (size(silabs, /type) eq 8) then begin optdep = spline(silabs.wavele, silabs.optdepth, wavele) - $ silabs.optdepth(0) absorp = exp(-optdep * sildep) endif else absorp = make_array(n_elements(nu), value=1.0) ; LOOP OVER WAVELENGTHS ;------------------------------------------------------------------------------- for i=0, n_elements(nu)-1 do begin ; CALCULATE FLUX OF DISK AT GIVEN POINT ;----------------------------------------------------------------------- flxval = 2.0 * h_const * nu[i]^3 / c_const^2 / $ (exp(h_const * nu[i] / k_const / temper) - 1) ; CONVERT RESULT TO JANSKYS AND ADD SILICATE ABSORPTION ;----------------------------------------------------------------------- flxval *= (2.0*!pi/(3600.*360.*1000./pscale))^2 * 1.0e26 * absorp[i] ; COPY RESULT INTO STRUCT ;----------------------------------------------------------------------- result[i].imgdat = flxval ; 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 ; PRINT TIME PROGRAMME TOOK TO RUN ;------------------------------------------------------------------------------- if not 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