;============================================================================ ;============================================================================ ; ; ZKERNEL ; ; DIRBE Kernel for use by ZMODEL to Fit the Zodiacal Dust Cloud ; ; Last Revision: 17 June 1996 (mark 6p8) ; ;============================================================================ ;============================================================================ ; ;---------------------------------------------------------------------------- ; Function GAUSSINT ; ; Returns Gauss-Laguerre abscissae and weights for 24 or 50 point gaussian ; quadrature. ; ; Written By: BA Franz, ARC, 9/94 ; Quadrature abscissae and weights computed by W. T. Reach ; ;---------------------------------------------------------------------------- pro gaussint,a,b,grid,wts,numpts=numpts if (n_elements(numpts) eq 0) then numpts = 24 ; ; Set-up static storage space common gaussint_cmn,init,grid24,wts24,grid25,wts25,grid50,wts50 ; ; Initialize abscissae and weights if first call ; if (n_elements(init) eq 0) then begin init = 1 grid24 = double( $ [0.0567047755,0.2990108986,0.7359095554,1.3691831160, $ 2.2013260537,3.2356758036,4.4764966151,5.9290837627, $ 7.5998993100,9.4967492209,11.629014912,14.007957977, $ 16.647125597,19.56289801, 22.775241987,26.308772391, $ 30.194291163,34.471097572,39.190608804,44.422349336, $ 50.264574993,56.864967174,64.466670616,73.534234792] $ )/73.534234792 wts24 = double( $ [0.1455497377,0.3393497722,0.5347365921,0.7322248725, $ 0.9326159014,1.1367925904,1.3457293379,1.560519046, $ 1.7824092263,2.0128491498,2.2535525026,2.5065825127, $ 2.7744704429,3.0603848697,3.3683805666,3.7037765832, $ 4.0737527888,4.4883345170,4.9621093140,5.5174318658, $ 6.19095,7.0486426246,8.2265276593,10.0758794961] $ )/73.534234792 grid50 = double( $ [0.0286305183,0.1508829357,0.3709487815,0.6890906999, $ 1.1056250235,1.6209617511,2.2356103759,2.9501833666, $ 3.7653997744,4.6820893876,5.7011975748,6.8237909098, $ 8.0510636694,9.3843453083,10.825109031,12.374981608, $ 14.035754599,15.809397197,17.698070933,19.704146535, $ 21.830223306,24.079151444,26.454057841,28.958376011, $ 31.595880956,34.370729963,37.287510610,40.351297573, $ 43.567720270,46.943043991,50.484267963,54.199244880, $ 58.096828017,62.187054175,66.481373878,70.992944826, $ 75.737011547,80.731404802,85.997211136,91.559690412, $ 97.449565614,103.70489123,110.37385880,117.51919820, $ 125.22547013,133.61202792,142.85832548,153.26037197, $ 165.3856433,180.69834370] $ )/180.69834370 wts50 = double( $ [0.0734786263,0.1711133062,0.2690583985,0.3672778840, $ 0.4658590232,0.5648992928,0.6644999757,0.7647657788, $ 0.8658052545,0.9677314398,1.0706625889,1.1747230029, $ 1.2800439495,1.3867646993,1.4950336890,1.6050098453, $ 1.7168640913,1.8307810796,1.9469611911,2.0656228595, $ 2.1870052872,2.3113716402,2.4390128272,2.5702520038, $ 2.7054499706,2.8450116969,2.9893942558,3.1391165624, $ 3.2947714220,3.4570405806,3.6267137126,3.8047126447, $ 3.9921226303,4.1902332693,4.4005928365,4.6250816069, $ 4.8660126478,5.1262732767,5.4095283377,5.7205203736, $ 6.0655271405,6.4530854743,6.8951889633,7.4093807809, $ 8.0226692310,8.7795282602,9.7603031266,11.131390723, $ 13.324267692,18.114146002] $ )/180.69834370 endif ; ; Select grid length based on user input ; case numpts of 50 : begin grid = grid50 wts = wts50 end else : begin grid = grid24 wts = wts24 endelse endcase ; ; Scale grid to user limits ; grid = temporary(grid) * (b-a) + a wts = temporary(wts) * (b-a) return end ; ;------------------------------------------------------------------------------- ; ; Procedure SimpInt ; ; Simpson Integration abscissae and weights as per HTF NEWESTZKERNEL0 ; Necessary for toroidal bands which are very narrow ; Added: JP 07 may 1996 ; pro simpint,a,b,grid,wts,stepsize=stepsize ; a is assumed to be zero! if (n_elements(stepsize) eq 0) then stepsize = 0.025 h = stepsize tsteps = round(b/h) + 1 grid = findgen(tsteps) * h wts = fltarr(tsteps) wts(0) = .333333 * h for i = 1,tsteps-1,2 do wts(i) = 1.333333 * h for i = 2,tsteps-2,2 do wts(i) = 0.666667 * h wts(tsteps-1) = .333333 * h return end ; ;------------------------------------------------------------------------------- ; ; Procedure EARTHSUN ; ; Returns Solar Elongation and Earth Position from Input Day and Lon, Lat ; ; Written By: BA Franz, WT Reach ; ; Inputs: ; Day - 1990 Day Number (1 = 01 Jan 1990) ; Lon - Ecliptic Lon of Line-of-Sight (Deg) ; Lat - Ecliptic Lat of Line-of-Sight (Deg) ;------------------------------------------------------------------------------- pro earthsun,day,lon,lat,SolElong,Earth_Dis,Earth_Lon,Earth_Mean_Lon pi2 = 2*!pi d2r = !pi/180. eccen = 0.01671254 lambda_solar = (-80.598349 + 0.98564736 * day + $ 1.912 * cos(pi2/365.25 * (day-94.8))) * d2r mean_anomaly = (356.637087 + 0.98560028 * day)*d2r mod pi2 Earth_Dis = (1.-eccen^2)/(1.+eccen*cos(mean_anomaly)) ; (AU) Earth_Lon = -(!pi-lambda_solar) mod (2.*!pi) ; (rad) SolElong = acos( cos(lat*d2r) * cos(lon*d2r - lambda_solar) ) Earth_Mean_Lon = ((99.403445 + 0.98564736*day)*d2r) mod pi2 return end ; ;------------------------------------------------------------------------------- ; ; Function COLCORR ; ; Returns DIRBE color correction coefficients for a modified BB source. ; ; Performs a table look-up of temperature versus correction. The color ; corrections are generated from the DIRBE Explanatory Supplement tables ; from Appendix B, for a black body with an emissivity index of 0. On first ; call, the Exp. Supp. data is interpolated to an even temperature grid ; and stored in an array as described below: ; ; Table(itemp,ival,idet) ; ; itemp = temperature index ; ival = value index (0=Color Correction, 1=dCdT) ; idet = detector index (0=1.25um, 9=240um) ; ; Written By: BA Franz, 10/95 ; ; Warning: Only applies to standard DIRBE wavelengths. ; ;------------------------------------------------------------------------------- function colcorr,det,t,dcdt ; ; Create static space for colcor correction tables and supporting data ; common colcor_cmn,init,temp,table,tmax,tmin,dt,nt ; ; If first call, initialize color correction tables ; if (n_elements(init) eq 0) then begin init = 1 ; ; Read Exp. Supp. Table for BB at P=0 ; ; Hardcode the file location ; findpro,'zkernel',/noprint,dirlist=dirlist ; if (dirlist(0) ne '') then begin ; case (!version.os) of ; 'vms' : infile = dirlist(0) + 'colcorr.tab' ; else : infile = dirlist(0) + '/' + 'colcorr.tab' ; endcase ; endif else infile = 'bafsci2:[dbswg.zl.zmodel.alpha]colcorr.tab' infile = 'colcorr.tab' print,'Reading color correction file: ',infile fmt = '(f,f,f,f,f,f,f,f,f,f,f)' readcol,infile,form=fmt,intemp,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10 incc = [[c1],[c2],[c3],[c4],[c5],[c6],[c7],[c8],[c9],[c10]] ; ; Define Interpolation Range ; tmax = 1000. tmin = 10. dt = 5. nt = (tmax - tmin)/dt + 1. temp = findgen(nt)*dt + tmin ; ; Build Table of CC and dCdT per Det ; table = fltarr(nt,2,10) for i=0,9 do begin cc = spline(intemp,incc(*,i),temp) cc2 = spline(intemp,incc(*,i),temp+dt) table(*,0,i) = cc table(*,1,i) = (cc2-cc)/dt endfor endif ; ; Compute CC and dCdT at Each Temperature in t ; it = fix((t - tmin)/dt) > 0 < (nt-1) delta_t = (t-temp(it)) dcdt = table(it,1,det) cc = table(it,0,det) + dcdt * delta_t return,cc end ; ;------------------------------------------------------------------------------- ;+ ; NAME: ; ZPLANCK ; ; PURPOSE: ; IDL function to return the spectral radiance of a blackbody, ; i.e. the Planck curve, in units of W/cm^2/sr. This is ; nu*I_nu (=lambda*I_lambda). The ; blackbody temperature and either frequency (in icm or GHz) ; or wavelength (in microns) are inputs to the function. The ; routine also optionally returns the derivative with respect to ; temperature, in units of W/cm^2/sr/K. ; ; CALLING SEQUENCE: ; result = planck (temperature, nu_or_lambda, [dBbT], [units=]) ; ; INPUTS: ; T float scalar Temperature of blackbody, in K. ; nu_or_lambda float [array] Frequency or wavelength at which to ; calculate spectrum. Units are as ; specified with "units" keyword. ; ; OPTIONAL INPUT: ; units keyword string 'Microns', 'icm', or 'GHz' to identify ; units of nu_or_lambda. Only first ; character is required. If left out, ; default is 'microns'. ; ; OUTPUTS: ; planck float [array] Spectral radiance at each wavelength. ; Units are W/cm^2/sr. ; ; OPTIONAL OUTPUT: ; dB/dT float [array] (Optional) derivative of Planck ; with respect to temperature. Units ; are W/cm^2/sr/K ; ; SIDE EFFECTS: ; None known. ; ; PROCEDURE: ; Identifies units using the "units" keyword, then converts the ; supplied independent variable into microns to evaluate the Planck ; function. Uses Rayleigh-Jeans and Wien approximations for the low- ; and high-frequency end, respectively. ; ; EXAMPLE: ; To produce a 35 K spectrum at 2,4,6,8,...,100 microns: ; ; wavelength = 2. + 2.*findgen(50) ; temp = 35. ; blackbody = planck (wavelength, temp, units='micron') ; ; One could also get back the derivative by including it in the call: ; blackbody = planck (wavelength, temp, deriv, units='m') ; ; RESTRICTIONS: ; Temperature must be given as a real*4 input, NOT an integer! ; Routine also gives incorrect results for T < 1 microkelvin (so sue me). ; ; REVISION HISTORY: ; Written by Rich Isaacman, General Sciences Corp. 17 April 1991 ; Allow array of Temperatures and single Lambda ; Rename to avoid confusion 6/95 ; ; REFERENCE: ; See Allen, Astrophysical Quantities for the Planck formula. ;- function zplanck, T, lambda, dBdT ; ; Define some necessary constants ; c = 2.9979e14 ;Speed of light (micron/sec) h = 6.6262e-34 ;Planck constant k = 1.3807e-23 ;Boltzmann constant stephan_boltzmann = 5.6697e-12 ;Stephan-Boltzmann constant coeff = 15. / !PI^5 * stephan_boltzmann ;Appropriate combination ; ; Introduce dimensionless variable chi, used to check whether we are on ; Wien or Rayleigh Jeans tails ; chi = h * c / (lambda * k * T) val = CHI*0 dBdT = CHI*0 ; Start on Rayleigh Jeans side ; rj = where (chi lt 0.001) if rj(0) ne -1 then begin val(rj) = coeff * T(rj)^4 * chi(rj)^3 dBdT(rj) = val(rj)/T(rj) endif ; ; Now do nonapproximate part ; exact = where (chi ge 0.001 and chi le 50) if exact(0) ne -1 then begin chi_ex = chi(exact) val(exact) = coeff * T(exact)^4 * chi_ex^4 / (exp(chi_ex) - 1.) dBdT(exact) = val(exact) * chi_ex / T(exact) / (1. - exp(-chi_ex)) endif ; ; ...and finally the Wien tail ; wien = where (chi gt 50.) if wien(0) ne -1 then begin chi_wn = chi(wien) val(wien) = coeff * T(wien)^4 * chi_wn^4 * exp(-chi_wn) dBdT(wien) = val(wien) * chi_wn / T(wien) / (1. - exp(-chi_wn)) endif ; if (n_elements(chi) eq 1) then val=val(0) if (n_elements(chi) eq 1) then dBdT=dBdT(0) return, val ; end ; ;-------------------------------------------------------------- ; Function THERMFUNC ; ; Returns the thermal component of the zodiacal dust source ; function. ; ; Written By: BA Franz, ARC, 9/95 ; ; Inputs: ; det - DIRBE detector number (0-9) or -1 if non-DIRBE ; Lambda - Wavelength in um ; R - sqrt(x^2+y^2+z^2) (AU) ; a - Model parameters ; ; Outputs: ; THERMFUNC - Brightness per particle in MJy/sr ; df - Partial derivative of source function wrt each ; parameter ; ;-------------------------------------------------------------- function thermfunc,det,Lambda,R,a,df,no_colcorr=no_colcorr if (n_params(0) gt 4) then want_partials=1 else want_partials=0 d2r = !pi/180. eps = double(1e-20) ; ; Thermal Parameters ; To1 = a(0) To2 = a(1) Kappa = a(2) Delta = a(3) ; ; Conversion from W/cm^2/Sr (Nu Bnu) to MJy/Sr (Bnu) ; Lambda/3e14 Hz^-1 * 1e4 cm^2/m^2 * 1e20 (MJy)/(W/m^2) ; cfact = Lambda / 3.0e-10 ; ; First Temperature Component ; Temp1 = To1/R^Delta Bnu1 = zPlanck(Temp1,Lambda,dBdT1) * cfact dBdT1 = dBdT1 * cfact if (det ge 0 and not keyword_set(no_colcorr)) then begin CCtherm1 = colcorr(det,temp1,dCdT1) endif else begin CCtherm1 = 1.0 dCdT1 = 0.0 endelse ; ; Second Temperature Component ; if (Kappa ne 0.0) then begin Temp2 = To2/R^Delta Bnu2 = zPlanck(Temp2,Lambda,dBdT2) * cfact dBdT2 = dBdT2 * cfact if (det ge 0 and not keyword_set(no_colcorr)) then begin CCtherm2 = colcorr(det,temp2,dCdT2) endif else begin CCtherm2 = 1.0 dCdT2 = 0.0 endelse endif else begin Temp2 = 0.0 Bnu2 = eps dBdT2 = 0.0 CCtherm2 = 1.0 dCdT2 = 0.0 endelse ; ; Total Thermal Contribution ; Therm = Bnu1*CCtherm1 + Kappa * Bnu2*CCtherm2 dBdT = dBdT1 + Kappa * dBdT2 ; ; Calculate Derivatives of Varying Parameters if Requested ; -------------------------------------------------------- if (want_partials) then begin npts = n_elements(R) npar = n_elements(a) if (n_elements(df) ne npts*npar) then df = dblarr(npts,npar) dT1 = dBdT1*CCTherm1 + Bnu1*dCdT1 dT2 = (dBdT2*CCTherm2 + Bnu2*dCdT2)*Kappa ; To1 df(*,0) = dT1/R^Delta ; To2 df(*,1) = dT2/R^Delta ; Kappa df(*,2) = Bnu2*CCtherm2 ; Delta df(*,3) = -1.*alog(R)*(Temp1*dT1 + Temp2*dT2) endif return,Therm end ; ;------------------------------------------------------------------------------- ; IDL Function PHASEFUNC.PRO ; ; Trial Phase Function for ZKERNEL ; ; Primary Author: JP, using BAF template Dec 1995 ; ; Contributors: ; ; Inputs: ; x - array of scattering angles (radians) ; a - Model parameters ; ; Optional Inputs: ; ; Outputs: ; phi - array of phase function ; pder - partial derivatives of phase function w.r.t each parameter ; ; Optional Outputs: ; ; Last modification: 18 March 1996 JP -- three parameter PF ; ;------------------------------------------------------------------------------- function phasefunc,x,a,pder C0 = a(0) C1 = a(1) ka = a(2) ; ; Evaluate the function ; pisq = !pi^2 picu = !pi^3 q0 = 2. q1 = !pi qk = (exp(ka*!pi) + 1.)/(ka*ka + 1) v = (2.*!pi)*( C0*q0 + C1*q1 + qk) u = C0 + C1*x + exp(ka*x) phi = u/v ; ; Evaluate the partial derivatives ; npts = n_elements(x) pder = dblarr(npts,3) pder(*,0) = (1.d0 - phi*4.*!pi)/v pder(*,1) = (x - 2.d0*!pi*phi*q1)/v pder(*,2) = (x*exp(ka*x) - 2.d0*!pi*phi*( (ka^2+1)*!pi*exp(ka*!pi) $ - 2.d0*ka*(exp(ka*!pi)+1))/((ka*ka+1.)^2))/v return,phi end ; ;-------------------------------------------------------------- ; Function SCATTFUNC ; ; Returns the scattered light component of the zodiacal dust ; source function. ; ; Written By: BA Franz, ARC, 11/95 ; Modified By: JP Dec 1995 ; ; Inputs: ; det - DIRBE detector number (0-9) or -1 if non-DIRBE ; LOS - Distance from Earth (AU) ; R - sqrt(x^2+y^2+z^2) (AU) ; Re - Distance from Earth to Sun (AU) ; SolElong - Solar elongation (radians) ; a - Model parameters ; ; Outputs: ; SCATTFUNC - Brightness per particle in MJy/sr ; df - Partial derivative of source function wrt each ; parameter ; ;-------------------------------------------------------------- function scattfunc,det,LOS,R,Re,SolElong,a,df if (n_params(0) gt 6) then want_partials=1 else want_partials=0 d2r = !pi/180. eps = double(1e-20) ; ; Solar Flux at 1 AU from Sun (courtesy W. T. Reach) ; SolFlux1AU = [ 2.3405606e+08, $ 1.2309874e+08, $ 64292872., $ 35733824., $ 5763843.0, $ 1327989.4, $ 230553.73, $ 82999.336, $ 42346.605, $ 14409.608 ] if (det ge 0) then begin SolFLux = SolFlux1AU(det) / R^2 endif else begin SolFlux = 0.0 endelse if (det le 2) then begin ; ; Evaluate Phase Function ; PhaseAng = asin( Re/R*sin(SolElong) < 1.0 ) itest = 1.0*(los ge Re*cos(SolElong)) ScatAng = (1.-itest)*PhaseAng + (itest)*(!pi - PhaseAng) aaa = a(3*det:3*det+2) Phi = phasefunc(ScatAng,aaa,dphi) ; ; Total Scattered Light Contribution ; CCscat = 1.D0 Scatt = SolFLux * CCscat * Phi endif else begin Scatt = 0.0d0 endelse ; ; Calculate Derivatives of Varying Parameters if Requested ; -------------------------------------------------------- if (want_partials) then begin npts = n_elements(R) npar = n_elements(a) if (n_elements(df) ne npts*npar) then df = dblarr(npts,npar) $ else df = temporary(df)*0.0 if (det le 2) then begin ; PF_C0 df(*,3*det) = SolFlux * CCscat * dphi(*,0) ; PF_C1 df(*,3*det+1) = SolFlux * CCscat * dphi(*,1) ; PF_Ka df(*,3*det+2) = SolFlux * CCscat * dphi(*,2) endif endif return,Scatt end ; ;-------------------------------------------------------------- ; Function ZSRCFUNC ; ; Returns the zodiacal dust source function, the brightness per ; dust grain. ; ; Written By: BA Franz, ARC, 9/95 ; ; Inputs: ; Scatt - Scattered light contribution ; Therm - Thermal emission contribution ; a - Model parameters (Emissivities/Albedos) ; ; Outputs: ; ZSRCFUNC - Brightness per particle in MJy/sr ; df - Partial derivative of source function wrt each ; parameter ; ;-------------------------------------------------------------- function zsrcfunc,det,Scatt,Therm,a,df,dScatt,dTherm if (n_params(0) gt 3) then want_partials=1 else want_partials=0 npar = n_elements(a) ; ; Albedo per Det ; AlbedoDet = [a(0:3),dblarr(6)+a(3)] ; ; Thermal Emissivity per Det ; EmissDet = [dblarr(2)+1.,a(4:11)] ; ; Set Albedo and Emissivity for This Detector ; ------------------------------------------- if (det ge 0) then begin Albedo = AlbedoDet(det) Emiss = EmissDet(det) endif else begin Albedo = 0. Emiss = 1. endelse ; ; Calculate Total Source per LOS Element ; ------------------------------------------ Source = Albedo * Scatt + Emiss * (1.0-Albedo) * Therm ; ; Calculate Derivatives of Varying Parameters if Requested ; -------------------------------------------------------- if (want_partials) then begin npts = n_elements(Therm) npar = n_elements(a) if (n_elements(df) ne npts*npar) then $ df = dblarr(npts,npar) $ else $ df = temporary(df)*0.0 dScatt = Albedo dTherm = Emiss * (1.0-Albedo) ; Albedo sign1 = 1.0*(Albedo ge 0) - (Albedo lt 0) sign2 = 1.0*((1.-Albedo) ge 0) - ((1.-Albedo) lt 0) dA = sign1*Scatt - sign2*Emiss*Therm df(*,(det < 3)) = dA ; Emissivity if (det ge 2) then $ df(*,4 + (det-2)) = Therm/Emiss endif return,Source end ; ;-------------------------------------------------------------- ; Function ZCloud ; ; Returns the number density of the zodiacal cloud. ; ; Written By: BA Franz, ARC, 9/95 ; ; Inputs: ; x,y,z - heliocentric cartesian coordinates ; R - sqrt(x^2+y^2+z^2) ; a - Model parameters ; ; Keyword Inputs: ; FuncIndex - selector for latitudinal density distribution ; 0 : Good ; 1 : Modified Fan ; 2 : Widened Modified Fan ; 3 : Ellipsoid ; 4 : Sombrero ; ; Outputs: ; ZCloud - number density ; df - partial of density wrt each parameter ; ;-------------------------------------------------------------- function zcloud,x,y,z,R,a,df,FuncIndx=FuncIndx if (n_params(0) gt 5) then want_partials=1 else want_partials=0 No = a(0) if (No ne 0.0) then begin d2r = !pi/180. Alpha = a(1) Beta = a(2) Gamma = a(3) Mu = a(4) Mu2 = a(5) Mu3 = a(6) Mu4 = a(7) Mu5 = a(8) Mu6 = a(9) Mu7 = a(10) Mu8 = a(11) Mu9 = a(12) Mu10 = a(13) Omega = a(14) * d2r Incl = a(15) * d2r Xo = a(16) Yo = a(17) Zo = a(18) ; ; Some frequently used geometry terms ; sino = sin(Omega) coso = cos(Omega) sini = sin(Incl) cosi = cos(Incl) ; ; Translate origin from Sun to cloud center. ; Xp = X - Xo Yp = Y - Yo Zp = Z - Zo Rc = sqrt( Xp^2 + Yp^2 + Zp^2 ) ; ; Rotate into coord system of cloud to determine Z-height in cloud ; (Zc) and radial distance in cloud symetry plane (Rc). ; Perform rotation of Omega around Z-Axis to align X with cloud ; line-of-nodes. Then, perform Incl rotation around X-Axis. ; Zc = sino*sini*Xp - coso*sini*Yp + cosi*Zp Zeta = abs(Zc/Rc) ; ; Compute Radial power-law parameter ; Mark6p2 sets this to Rc, not R 20 May 1996 ; Mark6p7 allows the polynomial mods only beyond a cutoff radius ; Mark6p8 puts is back the way it was, but keeps insistence on positive ; coeffs -- uses R, not Rc ; AlphaR = Alpha + ( Mu6^2*(R-1.) + Mu7^2*(R-1.)^2 + Mu8^2*(R-1.)^3 ) ;protection against overflow 12May1996 stuff = where(AlphaR gt 50.) if (stuff(0) ne -1) then AlphaR(stuff) = 50. ; ; Radial and Latitudinal Density Distribution ; case FuncIndx of ; ; John Good Density Function ; 0 : begin Dc = sqrt(Rc^2 - Zc^2) ZoD = abs(Zc/Dc) ZoD2G = ZoD^Gamma Dens_Cloud_Vert = exp( -Beta * ZoD2G ) Dens_Cloud_Rad = 1./Dc^AlphaR end ; ; Modified Fan Density Function ; 1 : begin Zeta2G = Zeta^Gamma Dens_Cloud_Vert = exp( -Beta * Zeta2G ) Dens_Cloud_Rad = 1./Rc^AlphaR end ; ; Widened Modified Fan Density Function ; 2 : begin GZR = (0.5*Zeta^2/Mu) * (Zeta lt Mu) + $ (Zeta-0.5*Mu) * (Zeta ge Mu) GZR2G = GZR^Gamma Dens_Cloud_Vert = exp( -Beta * GZR2G ) Dens_Cloud_Rad = 1./Rc^AlphaR end ; ; Ellipsoidal Density Function ; 3 : begin Bterm = (1. + (Beta*Zeta)^2) Dens_Cloud_Vert = 1./Bterm^Gamma Dens_Cloud_Rad = 1./Rc^AlphaR end ; ; Cosine (Sombrero) ; 4 : begin cosB = sqrt(1-Zeta^2) cosB2G = cosB^Gamma Dens_Cloud_Vert = (1. + Beta * cosB2G)/(1. + Beta) Dens_Cloud_Rad = 1./Rc^AlphaR end ; ; Modified Fan Density Function with Polynomial Z-height ; 5 : begin aZ = abs(Z) BetaZ = Beta + Mu3*aZ + Mu4*aZ^2 + Mu5*aZ^3 Dens_Cloud_Vert = exp( -BetaZ * Zeta ) Dens_Cloud_Rad = 1./Rc^AlphaR end endcase ; ; Total Density of Cloud Component ; Dens = No * Dens_Cloud_Vert * Dens_Cloud_Rad ; ; Partials of Density wrt Parameters ; if (want_partials) then begin npts = n_elements(x) npar = n_elements(a) if (n_elements(df) ne npts*npar) then df = dblarr(npts,npar) sZc = 1.0*(Zc ge 0.) - (Zc lt 0.) ; No df(*,0) = Dens/No case FuncIndx of 0 : lnR = alog(Dc) else : lnR = alog(Rc) endcase ; Alpha df(*,1) = -Dens*lnR ; Mu6 df(*,9) = -Dens*lnR*(R-1.) *2.*Mu6 ; Mu7 df(*,10) = -Dens*lnR*(R-1.)^2 *2.*Mu7 ; Mu8 df(*,11) = -Dens*lnR*(R-1.)^3 *2.*Mu8 case FuncIndx of ; ; John Good Density Function ; 0 : begin Zsqr = (1.-Zeta^2) dlnf = AlphaR*Zeta/Zsqr - $ Beta*Gamma*Zeta^(Gamma-1)*Zsqr^(-Gamma/2-1) ; Beta df(*,2) = -Dens*ZoD2G ; Gamma df(*,3) = -Beta*Dens*ZoD2G*alog(ZoD) end ; ; Modified Fan Density Function ; 1 : begin dlnf = -Beta*Gamma*Zeta^(Gamma-1) ; Beta df(*,2) = -Dens*Zeta2G ; Gamma df(*,3) = -Beta*Dens*Zeta2G*alog(Zeta) end ; ; Widened Modified Fan Density Function ; 2 : begin dlnf = -Beta*Gamma*GZR^(Gamma-1) * $ ( (Zeta/Mu) * (Zeta lt Mu) + 1.* (Zeta ge Mu) ) ; Beta df(*,2) = -Dens*GZR2G ; Gamma df(*,3) = -Beta*Dens*GZR2G*alog(GZR) ; Mu df(*,4) = Beta*Dens*0.5* Gamma * GZR^(Gamma-1)* $ ( (Zeta/Mu)^2 * (Zeta lt Mu) + (Zeta ge Mu) ) end ; ; Ellipsoidal Model ; 3 : begin dlnf = -2*Gamma*Beta^2*Zeta/(1+(Beta*Zeta)^2) ; Beta df(*,2) = -2.*Gamma/Beta*Dens*(1.-1./Bterm) ; Gamma df(*,3) = -Dens*alog(Bterm) end ; ; Cosine Model ; 4 : begin dlnf = -Gamma*Zeta*Beta/(1+Beta)/Dens_Cloud_Vert*cosB2G/cosB^2 ; Beta df(*,2) = Dens/Dens_Cloud_Vert* $ (-Dens_Cloud_Vert + cosB2G)/(1+Beta) ; Gamma df(*,3) = Dens/Dens_Cloud_Vert* $ (Beta*cosB2G/(1+Beta)*alog(cosB)) end ; ; Modified Fan Density Function ; 5 : begin dlnf = -BetaZ ; Beta df(*,2) = -Dens*Zeta ; Mu3 df(*,6) = -Dens*Zeta*aZ ; Mu4 df(*,7) = -Dens*Zeta*aZ^2 ; Mu5 df(*,8) = -Dens*Zeta*aZ^3 end endcase ; Omega df(*,14) = Dens*dlnf*sZc/Rc*(coso*sini*Xp+sino*sini*Yp)*d2r ; Incl df(*,15) = Dens*dlnf*sZc/Rc*(sino*cosi*Xp-coso*cosi*Yp-sini*Zp)*d2r ; Xo df(*,16) = Dens*(-dlnf*sZc/Rc*sini*sino + (dlnf*Zeta+AlphaR)*Xp/Rc^2) ; Yo df(*,17) = Dens*(dlnf*sZc/Rc*coso*sini + (dlnf*Zeta+AlphaR)*Yp/Rc^2) ; Zo df(*,18) = Dens*(-dlnf*sZc/Rc*cosi + (dlnf*Zeta+AlphaR)*Zp/Rc^2) endif endif else Dens = x*0 return,Dens end ; ;-------------------------------------------------------------- ; Function MigBand ; ; Returns the number density of a dust band. ; ; The model is a parameterization of the migrating band concept ; of W. T. Reach, as provided by Reach. ; ; Written By: BA Franz, ARC, 9/95 ; Modified By: JP 26 May 1996 based on HTF R cutoff idea ; changes meaning of dR, and makes p2r and Vr useless ; Modified By: JP 11 June 1996 to protect against floating underflow ; ; Inputs: ; x,y,z - heliocentric cartesian coordinates ; R - sqrt(x^2+y^2+z^2) ; a - Model parameters ; ; Outputs: ; migband - number density ; df - partial of density wrt each parameter ; ;-------------------------------------------------------------- function migband,x,y,z,R,a,df if (n_params(0) gt 5) then want_partials=1 else want_partials=0 No = a(0) if (No ne 0.0) then begin d2r = !pi/180. Dz = a(1) * d2r Dr = a(2) Ro = a(3) Vi = a(4) Vr = a(5) Pi = a(6) Pr = a(7) P2r = a(8) Omega = a(9) * d2r Incl = a(10) * d2r Xo = a(11) Yo = a(12) Zo = a(13) sino = sin(Omega) coso = cos(Omega) sini = sin(Incl) cosi = cos(Incl) ; ; Translate origin from Sun to cloud center. ; Xp = X - Xo Yp = Y - Yo Zp = Z - Zo ; ; Rotate into coord system of cloud to determine Z-height in cloud ; (Zc) and radial distance in cloud symetry plane (Rc). ; Perform rotation of Omega around Z-Axis to align X with cloud ; line-of-nodes. Then, perform Incl rotation around X-Axis. ; Zc = sino*sini*Xp - coso*sini*Yp + cosi*Zp Rc = sqrt( Xp^2 + Yp^2 + Zp^2 ) Zeta = abs(Zc/Rc) ZDz = Zeta / Dz RDr = Rc / Dr ViTerm = 1. + ZDz^Pi / Vi ; protect against underflow in exponential terms WtTerm = 1. + 0.*RDr Dens = 0.0 * Rc arg1 = RDr^20 arg2 = ZDz^6 ican = where ((arg1 le 86) and (arg2 le 86)) if (ican(0) ne -1) then begin WtTerm(ican) = 1. - exp(-arg1(ican)) Dens(ican) = No * (Ro/Rc(ican))^Pr * exp(-arg2(ican)) $ * Viterm(ican) * WtTerm(ican) endif ican = where((arg2 le 86) and (arg1 gt 86)) if (ican(0) ne -1) then $ Dens(ican) = No * (Ro/Rc(ican))^Pr * exp(-arg2(ican)) $ * Viterm(ican) ;WtTerm should = 1. ; if (want_partials) then begin npts = n_elements(x) npar = n_elements(a) if (n_elements(df) ne npts*npar) then df = dblarr(npts,npar) sZ = 1.0*(Zc ge 0.) - (Zc lt 0.) dlnf = (-6. * ZDz^5 + Pi/Vi * ZDz^(Pi-1) / ViTerm) / Dz ; No df(*,0) = Dens/No ; No ; Dz df(*,1) = -Dens*ZDz*dlnf*d2r ; Dr --- new ; now protect against underflow df(*,2) = 0. * Rc arg1 = ZDz^6 arg2 = RDr^20/Dr ican = where((arg1 le 86) and (arg2 le 86) and (Dens ne 0.0)) if (ican(0) ne -1) then $ df(ican,2) = No * (Ro/Rc(ican))^Pr * exp(-arg1(ican))*ViTerm(ican) $ * (-20 * RDr(ican)^20 * exp(-arg2(ican)) ) ; Ro ; Vi df(*,4) = -Dens*Vi^(-2) * ZDz^Pi / ViTerm ; Vr ; Pi df(*,6) = Dens * ZDz^Pi * alog(ZDz) / ViTerm / Vi ; Pr df(*,7) = Dens*alog(Ro/Rc) ; P2r ; Omega df(*,9) = Dens*dlnf*sZ/Rc*(coso*sini*Xp + sino*sini*Yp)*d2r ; Incl df(*,10) = Dens*dlnf*sZ/Rc*(sino*cosi*Xp - coso*cosi*Yp - sini*Zp)*d2r ; Xo df(*,11) = Dens*(-dlnf*sZ/Rc*sini*sino + (dlnf*Zeta+Pr)*Xp/Rc^2) ; Yo df(*,12) = Dens*(dlnf*sZ/Rc*coso*sini + (dlnf*Zeta+Pr)*Yp/Rc^2) ; Zo df(*,13) = Dens*(-dlnf*sZ/Rc*cosi + (dlnf*Zeta+Pr)*Zp/Rc^2) endif endif else Dens = x*0 return,Dens end ; ;-------------------------------------------------------------- ; Function SolRing ; ; Returns the number density of the Earth's resonant dust ring. ; The model is a gaussian toroidal ring with Earth-leading and ; Earth-trailing gaussian blobs. It is a parameterization of ; W. T. reach based on numerical images published by S. Dermott. ; ; Written By: BA Franz, ARC, 9/95 ; Modified By: JP, 01-Apr-1996 to allow 2nd radius thickness to ring ; 10-Apr-1996 to change R and Z laws in Dens_SR ; 24-May-1996 to change R law back to r^-2 from r-4 ; 11-Jun-1996 to avoid unecessary component calcs ; and protect against underflow ; 14-Jun-1996 speed up omega and incl derivs. ; ; ; ; Inputs: ; x,y,z - heliocentric cartesian coordinates ; R - sqrt(x^2+y^2+z^2) ; Theta - Earth mean longitude ; a - Model parameters ; ; Outputs: ; solring - number density ; df - partial of density wrt each parameter ; ;-------------------------------------------------------------- function solring,x,y,z,R,Theta,a,df if (n_params(0) gt 6) then want_partials=1 else want_partials=0 SR_No = a(0) LB_No = a(4) TB_No = a(10) if (SR_No ne 0.0 or LB_No ne 0.0 or TB_No ne 0) then begin d2r = !pi/180. ; ; Set Parameters ; ; Ring SR_R = a(1) SR_dR = a(2) ; SR_dR2 = a(21) ;JP -- outer thickness SR_dZ = a(3) ; Leading Blob LB_R = a(5) LB_dR = a(6) LB_Theta = a(7) * d2r LB_dTheta = a(8) * d2r LB_dZ = a(9) ; Trailing Blob TB_R = a(11) TB_dR = a(12) TB_Theta = a(13) * d2r TB_dTheta = a(14) * d2r TB_dZ = a(15) ; Symmetry plane Omega = a(16) * d2r Incl = a(17) * d2r Xo = a(18) Yo = a(19) Zo = a(20) ; ; Some frequently used geometry terms ; sino = sin(Omega) coso = cos(Omega) sini = sin(Incl) cosi = cos(Incl) ; ; Compute Density ; SR_Z = sino*sini*X - coso*sini*Y + cosi*Z dTheta = (atan(Y,X) - Theta) ; ; JP changes --no more use of SR_dR2 ; if (SR_No ne 0.0) then begin arg = (R-SR_R)^2/(SR_dR^2) + abs(SR_Z)/(SR_dZ) Dens_SR = 0. * R iloc = where (arg le 86) if (iloc(0) ne -1) then Dens_SR(iloc) = SR_No * exp(-arg(iloc)) endif else begin Dens_SR = 0.0 * R endelse ; ; JP mod - make blobs also exp(-z) ; if (LB_No ne 0.0) then begin LB_Delta = (dTheta-LB_Theta) mod (2*!pi) LB_Delta = LB_Delta + 2*!pi*((LB_Delta lt -!pi)-(LB_Delta gt !pi)) Dens_LB = 0. * R arg = (R-LB_R)^2/(LB_dR^2) $ +(LB_Delta)^2/(LB_dTheta^2) $ + abs(SR_Z)/(LB_dZ) iloc = where(arg le 86) if (iloc(0) ne -1) then Dens_LB(iloc) = LB_No * exp(-arg(iloc)) endif else begin Dens_LB = 0.0 * R endelse if (TB_No ne 0.0) then begin TB_Delta = (dTheta-TB_Theta) mod (2*!pi) TB_Delta = TB_Delta + 2*!pi*((TB_Delta lt -!pi)-(TB_Delta gt !pi)) Dens_TB = 0. * R arg =(R-TB_R)^2/(TB_dR^2) $ +(TB_Delta)^2/(TB_dTheta^2) $ +abs(SR_Z)/(TB_dZ) iloc = where(arg le 86) if (iloc(0) ne -1) then Dens_TB(iloc) = TB_No * exp(-arg(iloc)) endif else begin Dens_TB = 0.0 * R endelse ; ; end JP mod Dens = Dens_SR + Dens_LB + Dens_TB ; ; Partials of Density wrt Parameters ; if (want_partials) then begin npts = n_elements(x) npar = n_elements(a) if (n_elements(df) ne npts*npar) then df = dblarr(npts,npar) if (SR_No ne 0.0) then begin ; SR_No df(*,0) = Dens_SR/SR_No ; SR_R df(*,1) = Dens_SR * 2.*((R-SR_R))/(SR_dR^2) ;SR_dR df(*,2) = Dens_SR * 2.*((R-SR_R)^2)/(SR_dR^3) ; SR_dZ df(*,3) = Dens_SR * abs(SR_Z)/(SR_dZ^2) endif if (LB_No ne 0.0) then begin ; LB_No df(*,4) = Dens_LB/LB_No ; LB_R df(*,5) = Dens_LB * 2.*(R-LB_R)/LB_dR^2 ; LB_dR df(*,6) = Dens_LB * 2.*(R-LB_R)^2/LB_dR^3 ; LB_dTheta df(*,8) = Dens_LB * 2.*LB_Delta^2/LB_dTheta^3 * d2r ; LB_dZ df(*,9) = Dens_LB * abs(SR_Z)/LB_dZ^2 endif if (TB_No ne 0.0) then begin ; TB_No df(*,10) = Dens_TB/TB_No ; TB_R df(*,11) = Dens_TB * 2.*(R-TB_R)/TB_dR^2 ; TB_dR df(*,12) = Dens_TB * 2.*(R-TB_R)^2/TB_dR^3 ; TB_dTheta df(*,14) = Dens_TB * 2.*TB_Delta^2/TB_dTheta^3 * d2r ; TB_dZ df(*,15) = Dens_TB * abs(SR_Z)/TB_dZ^2 endif ; set up some common factors ; signZ = 1.0*(SR_Z ge 0.) - (SR_Z lt 0.) dfdz = -signZ * (Dens_SR/SR_dZ + Dens_LB/LB_dZ + Dens_TB/TB_dZ) ; RB_Omega df(*,16) = dfdz*(coso*sini*X + sino*sini*Y) * d2r ; RB_Incl df(*,17) = dfdz*(sino*cosi*X - coso*cosi*Y - sini*Z) * d2r endif endif else Dens = x*0 return,Dens end ; ;---------------------------------------------------------------------------- ; Function GET_PARNAMES ; ; Assign names to ZKERNEL parameters ; ; Written By: BA Franz, ARC, 2/94 ; Modified By : JP March 1996 to add parnames band 1,2,3 phase func ; JP April 1996 to add second ring radial law ; ;---------------------------------------------------------------------------- function get_parnames npars = 256 parname = strarr(npars) parname(*) = 'Unused' parname( 0) = 'FuncIndx' parname( 1) = 'PF1_C0' parname( 2) = 'PF1_C1' parname( 3) = 'PF1_Ka' parname( 4) = 'PF2_C0' parname( 5) = 'PF2_C1' parname( 6) = 'PF2_Ka' parname( 7) = 'PF3_C0' parname( 8) = 'PF3_C1' parname( 9) = 'PF3_Ka' parname( 10) = 'To1' parname( 11) = 'To2' parname( 12) = 'K' parname( 13) = 'Delta' parname( 14) = 'C_No' parname( 15) = 'C_Alpha' parname( 16) = 'C_Beta' parname( 17) = 'C_Gamma' parname( 18) = 'C_Mu' parname( 19) = 'C_Mu2' parname( 20) = 'C_Mu3' parname( 21) = 'C_Mu4' parname( 22) = 'C_Mu5' parname( 23) = 'C_Mu6' parname( 24) = 'C_Mu7' parname( 25) = 'C_Mu8' parname( 26) = 'C_Mu9' parname( 27) = 'C_Mu10' parname( 28) = 'C_Omega' parname( 29) = 'C_Incl' parname( 30) = 'C_Xo' parname( 31) = 'C_Yo' parname( 32) = 'C_Zo' parname( 33) = 'C_A1' parname( 34) = 'C_A2' parname( 35) = 'C_A3' parname( 36) = 'C_A4' parname( 37) = 'C_E3' parname( 38) = 'C_E4' parname( 39) = 'C_E5' parname( 40) = 'C_E6' parname( 41) = 'C_E7' parname( 42) = 'C_E8' parname( 43) = 'C_E9' parname( 44) = 'C_E10' parname( 45) = 'B1_No' parname( 46) = 'B1_Dz' parname( 47) = 'B1_Dr' parname( 48) = 'B1_Ro' parname( 49) = 'B1_Vi' parname( 50) = 'B1_Vr' parname( 51) = 'B1_Pi' parname( 52) = 'B1_Pr' parname( 53) = 'B1_P2r' parname( 54) = 'B1_Omega' parname( 55) = 'B1_Incl' parname( 56) = 'B1_Xo' parname( 57) = 'B1_Yo' parname( 58) = 'B1_Zo' parname( 59) = 'B1_A1' parname( 60) = 'B1_A2' parname( 61) = 'B1_A3' parname( 62) = 'B1_A4' parname( 63) = 'B1_E3' parname( 64) = 'B1_E4' parname( 65) = 'B1_E5' parname( 66) = 'B1_E6' parname( 67) = 'B1_E7' parname( 68) = 'B1_E8' parname( 69) = 'B1_E9' parname( 70) = 'B1_E10' parname( 71) = 'B2_No' parname( 72) = 'B2_Dz' parname( 73) = 'B2_Dr' parname( 74) = 'B2_Ro' parname( 75) = 'B2_Vi' parname( 76) = 'B2_Vr' parname( 77) = 'B2_Pi' parname( 78) = 'B2_Pr' parname( 79) = 'B2_P2r' parname( 80) = 'B2_Omega' parname( 81) = 'B2_Incl' parname( 82) = 'B2_Xo' parname( 83) = 'B2_Yo' parname( 84) = 'B2_Zo' parname( 85) = 'B2_A1' parname( 86) = 'B2_A2' parname( 87) = 'B2_A3' parname( 88) = 'B2_A4' parname( 89) = 'B2_E3' parname( 90) = 'B2_E4' parname( 91) = 'B2_E5' parname( 92) = 'B2_E6' parname( 93) = 'B2_E7' parname( 94) = 'B2_E8' parname( 95) = 'B2_E9' parname( 96) = 'B2_E10' parname( 97) = 'B3_No' parname( 98) = 'B3_Dz' parname( 99) = 'B3_Dr' parname(100) = 'B3_Ro' parname(101) = 'B3_Vi' parname(102) = 'B3_Vr' parname(103) = 'B3_Pi' parname(104) = 'B3_Pr' parname(105) = 'B3_P2r' parname(106) = 'B3_Omega' parname(107) = 'B3_Incl' parname(108) = 'B3_Xo' parname(109) = 'B3_Yo' parname(110) = 'B3_Zo' parname(111) = 'B3_A1' parname(112) = 'B3_A2' parname(113) = 'B3_A3' parname(114) = 'B3_A4' parname(115) = 'B3_E3' parname(116) = 'B3_E4' parname(117) = 'B3_E5' parname(118) = 'B3_E6' parname(119) = 'B3_E7' parname(120) = 'B3_E8' parname(121) = 'B3_E9' parname(122) = 'B3_E10' parname(123) = 'B4_No' parname(124) = 'B4_Dz' parname(125) = 'B4_Dr' parname(126) = 'B4_Ro' parname(127) = 'B4_Vi' parname(128) = 'B4_Vr' parname(129) = 'B4_Pi' parname(130) = 'B4_Pr' parname(131) = 'B4_P2r' parname(132) = 'B4_Omega' parname(133) = 'B4_Incl' parname(134) = 'B4_Xo' parname(135) = 'B4_Yo' parname(136) = 'B4_Zo' parname(137) = 'B4_A1' parname(138) = 'B4_A2' parname(139) = 'B4_A3' parname(140) = 'B4_A4' parname(141) = 'B4_E3' parname(142) = 'B4_E4' parname(143) = 'B4_E5' parname(144) = 'B4_E6' parname(145) = 'B4_E7' parname(146) = 'B4_E8' parname(147) = 'B4_E9' parname(148) = 'B4_E10' parname(149) = 'SR_No' parname(150) = 'SR_R' parname(151) = 'SR_dR' parname(152) = 'SR_dZ' parname(153) = 'LB_No' parname(154) = 'LB_R' parname(155) = 'LB_dR' parname(156) = 'LB_Theta' parname(157) = 'LB_dTheta' parname(158) = 'LB_dZ' parname(159) = 'TB_No' parname(160) = 'TB_R' parname(161) = 'TB_dR' parname(162) = 'TB_Theta' parname(163) = 'TB_dTheta' parname(164) = 'TB_dZ' parname(165) = 'RB_Omega' parname(166) = 'RB_Incl' parname(167) = 'RB_Xo' parname(168) = 'RB_Yo' parname(169) = 'RB_Zo' parname(170) = 'RB_A1' parname(171) = 'RB_A2' parname(172) = 'RB_A3' parname(173) = 'RB_A4' parname(174) = 'RB_E3' parname(175) = 'RB_E4' parname(176) = 'RB_E5' parname(177) = 'RB_E6' parname(178) = 'RB_E7' parname(179) = 'RB_E8' parname(180) = 'RB_E9' parname(181) = 'RB_E10' parname(182) = 'SR_dR2' ;JP return,parname end ; ;------------------------------------------------------------------------------- ; Function GET_PARINDEX ;------------------------------------------------------------------------------- function get_parindex,parname n = n_elements(parname) indx = intarr(n) namelist = strupcase(get_parnames()) for i=0,n-1 do begin name = strupcase(strtrim(parname(i),2)) ptr = where(namelist eq name) indx(i) = ptr(0) endfor if (n eq 1) then indx = indx(0) return,indx end ; ;------------------------------------------------------------------------------- ; Function WIRING2TOK ;------------------------------------------------------------------------------- function wiring2tok,wiring,input,output ptr1 = 0 ptr2 = strpos(wiring,'>',ptr1) if (ptr2 ne -1) then begin input = strtrim(strmid(wiring,ptr1,ptr2),2) len = strlen(wiring) nout = 0 done = 0 while (not done) do begin nout = nout+1 ptr1 = ptr2 + 1 ptr2 = strpos(wiring,'+',ptr1) if (ptr2 eq -1) then begin ptr2 = len done = 1 endif token = strtrim(strmid(wiring,ptr1,ptr2-ptr1),2) case nout of 1 : output = token else : output = [output,token] endcase endwhile endif return, nout end ; ;------------------------------------------------------------------------------- ; Function ZPAR_WIRING ; ; Builds the parameter wiring index from a list of wiring strings. If ; list is provided, it just returns the last installed index array. If ; no wiring was installed previously, and no wiring diagram was supplied, ; the function just returns a clean index array. ; ; BAF, 11/95 ;------------------------------------------------------------------------------- function zpar_wiring,wiring,reset=reset common zpar_wiring, wiring_indx npar = 256 ; ; If provided, convert wiring diagram to wiring index array ; if ( ((nset = n_elements(wiring))) gt 0 ) then begin ; ; Start with a cleared index ; wiring_indx = indgen(npar) ; ; Loop through each string, parsing to extract the input parameter ; index and the list of output parameter indices. ; for iset = 0,nset-1 do begin nout = wiring2tok(wiring(iset),input,output) if (nout ne 0) then begin ; Get indices indx_in = get_parindex(input) indx_out = get_parindex(output) ; Validate if (min([indx_in,indx_out]) lt 0) then begin message,'Invalid parameter name',/info stop endif ; Update the wiring index array. We must account for ; dependencies from previous strings, so we have to ; search the wiring index for each output index. for i=0,nout-1 do begin s = where( wiring_indx eq indx_out(i) ) wiring_indx(s) = indx_in endfor endif endfor endif else if (n_elements(wiring_indx) eq 0) then begin ; ; If it doesn't exist, create a dummy ; wiring_indx = indgen(npar) endif else if (keyword_set(reset)) then begin ; ; Clear the wiring index ; wiring_indx = indgen(npar) endif return,wiring_indx end ; ;------------------------------------------------------------------------------- ; Function PARFIX ; ; Function to fix parameters of ZKERNEL based on wavelength. ; ; Written By: BA Franz, ARC, 4/94 ; Modified By: JP March 1996 for wavelength-dependent phase function ; Modified By: JP May 26 1996 for migband changes ; ;------------------------------------------------------------------------------- function parfix,wl,a,indxpar npar = n_elements(a) indxa = indgen(npar) if (n_elements(indxpar) eq 0) then indxpar = indxa windx = zpar_wiring() ; ; Names of parameters which are currently inoperable or not optimizable ; pfix = 'FuncIndx' pfix = [pfix,'B1_Ro','B1_Vr','B1_P2r'] pfix = [pfix,'B2_Ro','B2_Vr','B2_P2r'] pfix = [pfix,'B3_Ro','B3_Vr','B3_P2r'] pfix = [pfix,'B4_Ro','B4_Vr','B4_P2r'] pfix = [pfix,'LB_Theta'] pfix = [pfix,'TB_Theta'] pfix = [pfix,'RB_Xo','RB_Yo','RB_Zo'] ; ; Names of parameters which require a specific wavelength ; s = where(wl eq 1.25) & if (s(0) eq -1) then begin pfix = [pfix,'C_A1','B1_A1','B2_A1','B3_A1','B4_A1','RB_A1'] pfix = [pfix,'PF1_C0','PF1_C1','PF1_Ka'] endif s = where(wl eq 2.20) & if (s(0) eq -1) then begin pfix = [pfix,'C_A2','B1_A2','B2_A2','B3_A2','B4_A2','RB_A2'] pfix = [pfix,'PF2_C0','PF2_C1','PF2_Ka'] endif s = where(wl eq 3.50) & if (s(0) eq -1) then begin pfix = [pfix,'C_A3','B1_A3','B2_A3','B3_A3','B4_A3','RB_A3'] pfix = [pfix,'C_E3','B1_E3','B2_E3','B3_E3','B4_E3','RB_E3'] pfix = [pfix,'PF3_C0','PF3_C1','PF3_Ka'] endif s = where(wl eq 4.90) & if (s(0) eq -1) then begin pfix = [pfix,'C_A4','B1_A4','B2_A4','B3_A4','B4_A4','RB_A4'] pfix = [pfix,'C_E4','B1_E4','B2_E4','B3_E4','B4_E4','RB_E4'] endif s = where(wl eq 12.) & if (s(0) eq -1) then begin pfix = [pfix,'C_E5','B1_E5','B2_E5','B3_E5','B4_E5','RB_E5'] endif s = where(wl eq 25.) & if (s(0) eq -1) then begin pfix = [pfix,'C_E6','B1_E6','B2_E6','B3_E6','B4_E6','RB_E6'] endif s = where(wl eq 60.) & if (s(0) eq -1) then begin pfix = [pfix,'C_E7','B1_E7','B2_E7','B3_E7','B4_E7','RB_E7'] endif s = where(wl eq 100.) & if (s(0) eq -1) then begin pfix = [pfix,'C_E8','B1_E8','B2_E8','B3_E8','B4_E8','RB_E8'] endif s = where(wl eq 140.) & if (s(0) eq -1) then begin pfix = [pfix,'C_E9','B1_E9','B2_E9','B3_E9','B4_E9','RB_E9'] endif s = where(wl eq 240.) & if (s(0) eq -1) then begin pfix = [pfix,'C_E10','B1_E10','B2_E10','B3_E10','B4_E10','RB_E10'] endif ; ; Convert names to indices of parameters we want fixed ; pfixIndx = get_parindex(pfix) ; ; Fix params which are unused ; s = where(a eq -1) if (s(0) ne -1) then pfixIndx = [pfixIndx,s] ; ; Fix params for components which are turned-off ; s = get_parindex(['K','To2']) if (a(s(0)) eq 0.0) then $ pfixIndx = [pfixIndx,s] s = get_parindex('B1_No') if (a(s(0)) eq 0.0) then $ pfixIndx = [pfixIndx,s(0)+indgen(26)] s = get_parindex('B2_No') if (a(s(0)) eq 0.0) then $ pfixIndx = [pfixIndx,s(0)+indgen(26)] s = get_parindex('B3_No') if (a(s(0)) eq 0.0) then $ pfixIndx = [pfixIndx,s(0)+indgen(26)] s = get_parindex('B4_No') if (a(s(0)) eq 0.0) then $ pfixIndx = [pfixIndx,s(0)+indgen(26)] s = get_parindex('SR_No') if (a(s(0)) eq 0.0) then $ pfixIndx = [pfixIndx,s(0)+indgen(4)] s = get_parindex('LB_No') if (a(s(0)) eq 0.0) then $ pfixIndx = [pfixIndx,s(0)+indgen(6)] s = get_parindex('TB_No') if (a(s(0)) eq 0.0) then $ pfixIndx = [pfixIndx,s(0)+indgen(6)] s = get_parindex(['SR_No','LB_No','TB_No']) if (max(a(s)) eq 0.0) then $ pfixIndx = [pfixIndx,s(0)+indgen(33)] ; ; Enforce unique, sorted order ; if (n_elements(pfixIndx) gt 1) then begin pfixIndx = pfixIndx(uniq(pfixIndx,sort(pfixIndx))) endif if (n_elements(indxpar) gt 1) then begin indxpar = indxpar(uniq(indxpar,sort(indxpar))) endif ; ; Apply pfixIndx to indxpar ; s = nowhere(indxpar,pfixIndx) if (s(0) ne -1) then indxpar = indxpar(s) ; ; Fix parameters which are wired to other parameters ; s = anywhere(indxpar,windx) if (s(0) ne -1) then $ indxpar = indxpar(s) $ else $ print,'Warning: parfix found no free parameters' return,indxpar end ; ;------------------------------------------------------------------------------- ; IDL Function ZKERNEL.PRO ; ; Zodiacal Dust Model ; ; Primary Authors: B.A. Franz, ARC, (301) 513-7776 ; W.T. Reach, USRA (301) 286-4255 ; Modified By: JP March 1996 for phase function/scattering function ; JP April 1996 for second radial law ; JP April 1996 as per HTF 4/16/96 code to include max. dist ; from sun for integration limits ; HTF CODE IN CAPS. ; ; Inputs: ; x - n-element structured array of independent variables which at least ; contains the following fields. (See ZODI_STR.PRO) ; x.longitude : Heliocentric Ecliptic longitude of LOS in degrees ; x.latitude : Heliocentric Ecliptic latitude of LOS in degrees ; x.day1990 : 1990 Day Number of Observation ; x.wave_len : Wavelength in um ; ; a - m-element floating array of model parameters ; ; Optional Inputs: ; indxpar=indxpar - array of indices of variable parameters (Def=0..m-1) ; ; Outputs: ; f - n-element floating array of Model Flux in MJy/Sr ; ; Optional Outputs: ; df - nxm array of partial derivatives for model parameters at each x ; ;------------------------------------------------------------------------------- ; pro zkernel,data,a,f,df,indxpar=indxpar,losinfo=losinfo,no_colcorr=no_colcorr if (keyword_set(losinfo)) then want_los_info=1 else want_los_info=0 ; ; Define some useful params ; RMAX = 5.2 ; JUPITER'S ORBIT RADIUS NUMBER_OF_STEPS = 50 ; THIS MANY FUNCTION EVALUATIONS PER L.O.S. npts = n_elements(data) ; d2r = !pi/180.D ; degrees to radians eps = double(1e-20) ; > zero dbwave = [1.25,2.2,3.5,4.9,12.,25.,60.,100.,140.,240.] nmsg = 50000 ; ; Set wiring ; windx = zpar_wiring() a = a(windx) ; ; Set density function selection flag ; FuncIndx = fix(a(0)) ; ; Get Scatt Function Parameters (JP) ; nScatt = 9 iScatt = 1 + indgen(nScatt) aScatt = a(iScatt) ; ; Get Therm Function Parameters ; nTherm = 4 iTherm = 10 + indgen(nTherm) aTherm = a(iTherm) ; ; Get Cloud Parameters ; nDens_C = 19 ; Number of density params iDens_C = 14 + indgen(nDens_C) ; Index of density params in a aDens_C = a(iDens_C) ; Vector of density params nSrc_C = 12 ; Number of emissivity params iSrc_C = 33 + indgen(nSrc_C) ; Index of emissivity params in a aSrc_C = a(iSrc_C) ; Vector of emissivity params ; ; Get Dust Band 1 Parameters ; nDens_B1 = 14 ; Number of density params iDens_B1 = 45 + indgen(nDens_B1) ; Index of density params in a aDens_B1 = a(iDens_B1) ; Vector of density params nSrc_B1 = 12 ; Number of emissivity params iSrc_B1 = 59 + indgen(nSrc_B1) ; Index of emissivity params in a aSrc_B1 = a(iSrc_B1) ; Vector of emissivity params ; ; Get Dust Band 2 Parameters ; nDens_B2 = 14 ; Number of density params iDens_B2 = 71 + indgen(nDens_B2) ; Index of density params in a aDens_B2 = a(iDens_B2) ; Vector of density params nSrc_B2 = 12 ; Number of emissivity params iSrc_B2 = 85 + indgen(nSrc_B2) ; Index of emissivity params in a aSrc_B2 = a(iSrc_B2) ; Vector of emissivity params ; ; Get Dust Band 3 Parameters ; nDens_B3 = 14 ; Number of density params iDens_B3 = 97 + indgen(nDens_B3) ; Index of density params in a aDens_B3 = a(iDens_B3) ; Vector of density params nSrc_B3 = 12 ; Number of emissivity params iSrc_B3 = 111 + indgen(nSrc_B3) ; Index of emissivity params in a aSrc_B3 = a(iSrc_B3) ; Vector of emissivity params ; ; Get Dust Band 4 Parameters ; nDens_B4 = 14 ; Number of density params iDens_B4 = 123 + indgen(nDens_B4) ; Index of density params in a aDens_B4 = a(iDens_B4) ; Vector of density params nSrc_B4 = 12 ; Number of emissivity params iSrc_B4 = 137 + indgen(nSrc_B4) ; Index of emissivity params in a aSrc_B4 = a(iSrc_B4) ; Vector of emissivity params ; ; Get Solar Ring Parameters ; nDens_RB = 21 ; Number of density params iDens_RB = 149 + indgen(nDens_RB) ; Index of density params in a nDens_RB = nDens_RB + 1 ; JP add in SR_dR2 iDens_RB = [iDens_RB,182] ; JP add in SR_dR2 aDens_RB = a(iDens_RB) ; Vector of density params nSrc_RB = 12 ; Number of emissivity params iSrc_RB = 170 + indgen(nSrc_RB) ; Index of emissivity params in a aSrc_RB = a(iSrc_RB) ; Vector of emissivity params ; ; Associate wavelengths with fullband detector numbers ; 1a=1b=1c=0 .... 9=b10 ; detnum = bytarr(npts) - 1 for i=0,9 do $ detnum = temporary(detnum) + ( data.wave_len eq dbwave(i) ) * (i+1) ; ; Calculate position of Earth in heliocentric ecliptic coords and ; Solar elongation of LOS ; earthsun,data.day1990,data.longitude,data.latitude, $ SolElong,Earth_Dis,Earth_Lon,Earth_Mean_Lon ; THIS PART GETS COMMENTED OUT WHEN CHANGE TO HELIOCENTRIC INTEGRATION RANGES ;------------------------------------------------ ; Set grid, weights for Gauss-Leguerre quadrature ; Define LOS Integration Points (from Earth in AU) ; ;gaussint,0.,5.,los,gqwts,numpts=50 ;------------------------------------------------ ; ; Set-up for calculation of integrated flux ; if (n_elements(f) ne npts) then f = dblarr(npts) ; ; Set-up for calculation of partial derivatives ; if (n_params(0) gt 3) then begin want_partials = 1 ; ; Create index of variable parameters, if not supplied. ; npar = n_elements(a) if (n_elements(indxpar) eq 0) then begin nterms = npar indxpar = indgen(nterms) endif else $ nterms = n_elements(indxpar) ; ; Create storage for partials, and create mask to expedite selection ; of parameters which require partials. ; if (n_elements(df) ne npts*nterms) then $ df = dblarr(npts,nterms) $ else $ df = df * 0.0 parnum = intarr(npar)-1 & parnum(indxpar) = indgen(nterms) endif else want_partials = 0 ; ; Set-up to store line-of-sight info ; if (want_los_info) then begin element_rec = { elemnt_str, $ earth_dist : 0.0, $ solar_dist : 0.0, $ ecliptic_z : 0.0, $ density : 0.0, $ flux : 0.0, $ int_wts : 0.0 $ } losdata = replicate( $ { wave_len : 0.0, $ pixel_no : 0L, $ lon : 0.0, $ lat : 0.0, $ elong : 0.0, $ day : 0.0, $ earth_lon : 0.0, $ earth_rad : 0.0, $ element : replicate( {elemnt_str},n_elements(los) ), $ zodi : 0.0 }, npts) endif ; ; Loop over each LOS, Time, Wavelength ; for ilos = 0L,npts-1 do begin ; ; Compute Basic Geometry ; ---------------------- lat = double(data(ilos).latitude * d2r ) lon = double(data(ilos).longitude * d2r ) ; ; Position of Earth ; Re = Earth_Dis(ilos) Theta = Earth_Lon(ilos) COSTHETA = COS(THETA) SINTHETA = SIN(THETA) COSLON = COS(LON) SINLON = SIN(LON) COSLAT = COS(LAT) SINLAT = SIN(LAT) X0 = RE*COSTHETA & Y0 = RE*SINTHETA ; &Z0 = 0.0 B = 2.*( X0*COSLAT*COSLON + Y0*COSLAT*SINLON ) C = RE*RE - RMAX*RMAX Q = -0.5*B*(1.0 + SQRT( B*B - 4.*C)/ABS(B)) RANGE = MAX( [Q, C/Q] ) ; Set up Integration abscissae and wts, depending upon ecliptic latitude ; Change by JP 08 May 1996 set to 20 rather than 15 ; betathresh = 20. * d2r ; ; THIS PART GETS INSIDE THE LOOP WHEN CHANGE TO HELIOCENTRIC INTEGRATION RANGES ;------------------------------------------------ ; Set grid, weights for Gauss-Leguerre quadrature ; Define LOS Integration Points (from Earth in AU) ; if (abs(lat) gt betathresh) then begin gaussint,0.,RANGE,los,gqwts,numpts=NUMBER_OF_STEPS endif else begin simpint,0., RANGE,los,gqwts,stepsize=0.025 endelse ;------------------------------------------------ ; ; Position of LOS elements in Heliocentric Ecliptic coord system. ; Sxy = los*COSLAT X = Re*COSTHETA + Sxy*COSLON Y = Re*SINTHETA + Sxy*SINLON Z = los*SINLAT R = sqrt( X^2 + Y^2 + Z^2 ) ; ; Source Function per LOS element ; ------------------------------------------------------------------ Lambda = double(data(ilos).wave_len) Det = detnum(ilos) Scatt = scattfunc(Det,los,R,Re,SolElong(ilos),aScatt,dScatt) Therm = thermfunc(Det,Lambda,R,aTherm,dTherm,no_colcorr=no_colcorr) ; ; Calculate Density and Source Associated With Smooth Cloud Component ; ------------------------------------------------------------------- Dens_C = zcloud(x,y,z,R,aDens_C,dDens_C,func=FuncIndx) Src_C = zsrcfunc(Det,Scatt,Therm,aSrc_C,dSrc_C,dSrc_dScatt_C,dSrc_dTherm_C) ; ; Calculate Density and Source Associated with Dust Bands ; ------------------------------------------------------------------- Dens_B1 = migband(x,y,z,R,aDens_B1,dDens_B1) Src_B1 = zsrcfunc(Det,Scatt,Therm,aSrc_B1,dSrc_B1,dSrc_dScatt_B1,dSrc_dTherm_B1) Dens_B2 = migband(x,y,z,R,aDens_B2,dDens_B2) Src_B2 = zsrcfunc(Det,Scatt,Therm,aSrc_B2,dSrc_B2,dSrc_dScatt_B2,dSrc_dTherm_B2) Dens_B3 = migband(x,y,z,R,aDens_B3,dDens_B3) Src_B3 = zsrcfunc(Det,Scatt,Therm,aSrc_B3,dSrc_B3,dSrc_dScatt_B3,dSrc_dTherm_B3) Dens_B4 = migband(x,y,z,R,aDens_B4,dDens_B4) Src_B4 = zsrcfunc(Det,Scatt,Therm,aSrc_B4,dSrc_B4,dSrc_dScatt_B4,dSrc_dTherm_B4) ; ; Calculate Density and Source Associated with Solar Ring/Blobs ; ------------------------------------------------------------------- Dens_RB = solring(x,y,z,R,Earth_Mean_Lon(ilos),aDens_RB,dDens_RB) Src_RB = zsrcfunc(Det,Scatt,Therm,aSrc_RB,dSrc_RB,dSrc_dScatt_RB,dSrc_dTherm_RB) ; ; Calculate Total Zodi Brightness per LOS Element ; -------------------------------------------------------------------- Flux = Src_C * Dens_C + $ Src_B1 * Dens_B1 + $ Src_B2 * Dens_B2 + $ Src_B3 * Dens_B3 + $ Src_B4 * Dens_B4 + $ Src_RB * Dens_RB ; ; Integrate Along LOS to Get Total Intensity ; ------------------------------------------ f(ilos) = total(gqwts*Flux) ; ; Store line-of-sight info (Needed by ZODI Explorer) ; -------------------------------------------------- if (want_los_info) then begin Dens = Dens_C + Dens_B1 + Dens_B2 + Dens_B3 + Dens_RB losdata(ilos).wave_len = Lambda losdata(ilos).pixel_no = data(ilos).pixel_no losdata(ilos).lon = Lon / d2r losdata(ilos).lat = Lat / d2r losdata(ilos).elong = SolElong(ilos) / d2r losdata(ilos).day = data(ilos).day1990 losdata(ilos).earth_lon = Theta / d2r losdata(ilos).earth_rad = Re losdata(ilos).element.earth_dist = LOS losdata(ilos).element.solar_dist = R losdata(ilos).element.ecliptic_z = Z losdata(ilos).element.density = Dens losdata(ilos).element.flux = Flux losdata(ilos).element.int_wts = gqwts losdata(ilos).zodi = f(ilos) endif ; ; Calculate Derivatives of Varying Parameters if Requested ; -------------------------------------------------------- if (want_partials) then begin ; ; Scatt Function Partials ; for ii = 0,nScatt-1 do begin ipar = parnum(windx(iScatt(ii))) if (ipar ge 0) then begin dSrcS = dScatt(*,ii) * ( $ dSrc_dScatt_C * Dens_C + $ dSrc_dScatt_B1 * Dens_B1 + $ dSrc_dScatt_B2 * Dens_B2 + $ dSrc_dScatt_B3 * Dens_B3 + $ dSrc_dScatt_B4 * Dens_B4 + $ dSrc_dScatt_RB * Dens_RB ) df(ilos,ipar) = df(ilos,ipar) + total(gqwts*dSrcS) endif endfor ; ; Therm Function Partials ; for ii = 0,nTherm-1 do begin ipar = parnum(windx(iTherm(ii))) if (ipar ge 0) then begin dSrcT = dTherm(*,ii) * ( $ dSrc_dTherm_C * Dens_C + $ dSrc_dTherm_B1 * Dens_B1 + $ dSrc_dTherm_B2 * Dens_B2 + $ dSrc_dTherm_B3 * Dens_B3 + $ dSrc_dTherm_B4 * Dens_B4 + $ dSrc_dTherm_RB * Dens_RB ) df(ilos,ipar) = df(ilos,ipar) + total(gqwts*dSrcT) endif endfor ; ; Cloud Partials ; if (aDens_C(0) ne 0.0) then begin for ii = 0,nDens_C-1 do begin ipar = parnum(windx(iDens_C(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*Src_C*dDens_C(*,ii)) endfor for ii = 0,nSrc_C-1 do begin ipar = parnum(windx(iSrc_C(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*dSrc_C(*,ii)*Dens_C) endfor endif ; ; Band 1 Partials ; if (aDens_B1(0) ne 0.0) then begin for ii = 0,nDens_B1-1 do begin ipar = parnum(windx(iDens_B1(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*Src_B1*dDens_B1(*,ii)) endfor for ii = 0,nSrc_B1-1 do begin ipar = parnum(windx(iSrc_B1(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*dSrc_B1(*,ii)*Dens_B1) endfor endif ; ; Band 2 Partials ; if (aDens_B2(0) ne 0.0) then begin for ii = 0,nDens_B2-1 do begin ipar = parnum(windx(iDens_B2(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*Src_B2*dDens_B2(*,ii)) endfor for ii = 0,nSrc_B2-1 do begin ipar = parnum(windx(iSrc_B2(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*dSrc_B2(*,ii)*Dens_B2) endfor endif ; ; Band 3 Partials ; if (aDens_B3(0) ne 0.0) then begin for ii = 0,nDens_B3-1 do begin ipar = parnum(windx(iDens_B3(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*Src_B3*dDens_B3(*,ii)) endfor for ii = 0,nSrc_B3-1 do begin ipar = parnum(windx(iSrc_B3(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*dSrc_B3(*,ii)*Dens_B3) endfor endif ; ; Band 4 Partials ; if (aDens_B4(0) ne 0.0) then begin for ii = 0,nDens_B4-1 do begin ipar = parnum(windx(iDens_B4(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*Src_B4*dDens_B4(*,ii)) endfor for ii = 0,nSrc_B4-1 do begin ipar = parnum(windx(iSrc_B4(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*dSrc_B4(*,ii)*Dens_B4) endfor endif ; ; Solar Ring/Blob Partials ; if (aDens_RB(0) ne 0.0) then begin for ii = 0,nDens_RB-1 do begin ipar = parnum(windx(iDens_RB(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*Src_RB*dDens_RB(*,ii)) endfor for ii = 0,nSrc_RB-1 do begin ipar = parnum(windx(iSrc_RB(ii))) if (ipar ge 0) then $ df(ilos,ipar) = df(ilos,ipar) $ + total(gqwts*dSrc_RB(*,ii)*Dens_RB) endfor endif endif if ((ilos+1) mod nmsg eq 0) then print,'LOS # ',ilos endfor ; Loop over LOS ; ; Replace returned flux with LOS structure if requested ; if (want_los_info) then f = temporary(losdata) return end