;+ ; NAME: ; rt_random ; PURPOSE: (one line) ; calculate molecular opacity in a specified frequency range ; DESCRIPTION: ; ********************************************************************* ; * * ; * RANDOM calculate the molecular opacity in a specified * ; * frequency range. Required input consists of an array of * ; * line positions, widths, strengths, and quantum numbers. * ; * The S-1 random model with a Curtis-Godson approximation * ; * is used. The relevant equations are described by Wallace, * ; * Prather, and Belton, Ap. J., 1974, Vol 193, pp 481-493. * ; * * ; ********************************************************************* ; ; CATEGORY: ; RT ; CALLING SEQUENCE: ; a = rt_random (uu,prs,tmp,col,str,wid,eng,efc,rmass,wmid,delnu) ; INPUTS: ; INPUT ; ; uu - cosine of solar zenith angle ; prs - pressure grid (atm) ; not ubar? [LAY] ; tmp - temperature grid (Kelvins) ; col - column density grid (cm-2) ; str - list of line strengths (cm-1) ; wid - Lorentz widths of lines at a ref. temperature ; (cm-1) ; and std pressure? [LAY] ; eng - energy of lower state (cm-1) ; rmass - molecular mass in g ; efc - rotational partition function exponent ; wmid - frequency in center of interval (cm-1) ; delnu - range of frequency covered by these lines ; OPTIONAL INPUT PARAMETERS: ; none ; KEYWORD INPUT PARAMETERS: ; none ; KEYWORD OUTPUT PARAMETERS: ; none ; OUTPUTS: ; a - average absorption in this frequency interval ; COMMON BLOCKS: ; SIDE EFFECTS: ; RESTRICTIONS: ; None ; PROCEDURE: ; ; ; MODIFICATION HISTORY: ; Written 2004 May 12, Leslie Young SwRI ; Modified from random.f, Roger Yelle ;- function rt_random, uu,prs,tmp,col,str,wid,eng,efc,rmass,wmid,delnu ntau = n_elements(prs) nlines = n_elements(str) hck = !phys.h * !phys.c/!phys.k if (wmid gt 2809.030 and wmid lt 2818.894) and $ (uu gt 0.1 and uu lt 0.12) then begin idebug=1 print, 'DEBUG random' endif else begin idebug=0 end if idebug eq 1 then begin print, 'str' print, str[0:2], form='(3e11.3)' print, '..', str[nlines-2:nlines-1], form='(A11,2e11.3)' print, 'wid' print, wid[0:2], form='(3e11.3)' print, '..', wid[nlines-2:nlines-1], form='(A11,2e11.3)' print, 'eng' print, eng[0:2], form='(3e11.3)' print, '..', str[nlines-2:nlines-1], form='(A11,2e11.3)' end ; ------------------------------------------------------------------ ; Internal arrays ; ------------------------------------------------------------------ ; dimension wdop(nd1),wlor(nd1),hckt(nd1),qrat(nd1),sumstr(nd1) ; & ,sumsqr(nd1),sumlor(nd1),u(nd1),y(nd1),f(nd1),alor(nd1),w(nd1) ; & ,v(nd1),g(nd1),adop(nd1),st(nd1,500),dcol(nd1) tmp0 = 296.d0 twopi = 2.d*!pi a = dblarr(ntau) if (nlines eq 0) then begin return, a endif wdop = wmid*sqrt(2.*!phys.k*tmp/rmass)/!phys.c wlor = 1.d-6*prs * sqrt(tmp0 / tmp) hckt = hck*(1.d/tmp - 1.d/tmp0) qrat = ((tmp0/tmp)^efc) ; -------------------------------------------------------------- ; Employ the Curtis-Godson approximation to determine ; appropriate averages for the pressure and column ; abundance. Use the local temperature as the reference. ; -------------------------------------------------------------- st = dblarr(ntau, nlines) for lw = 0, nlines-1 do begin for n = 0, ntau-1 do begin st[n,lw]=str[lw]*qrat[n]*exp(-hckt[n]*eng[lw]) endfor endfor dcol = dblarr(ntau) for n = 1, ntau-1 do begin dcol[n] = col[n]-col[n-1] endfor sumstr = dblarr(ntau) sumsqr = dblarr(ntau) sumlor = dblarr(ntau) for n = 0, ntau-1 do begin ; ----------------------------------------------------------- ; Calculate sum of line strengths ; ----------------------------------------------------------- suml = 0.d sumw = 0.d for lw = 0, nlines-1 do begin s1=st[0,lw] sum = s1*col[0] for m = 1, n do begin s2=st[m,lw] sum = sum + 0.5d*(s1+s2)*dcol[m] s1 = s2 endfor suml = suml + sum/uu sumw = sumw + sqrt(sum/uu) endfor sumstr[n] = suml sumsqr[n] = sumw ; ----------------------------------------------------------- ; Calculate sum of line strengths weighted by Lorentz width ; ----------------------------------------------------------- suml = 0. for lw = 0, nlines-1 do begin s1=st[0,lw]*wid[lw]*wlor[0] sum = s1*col[0] for m = 1, n do begin s2=st[m,lw]*wid[lw]*wlor[m] sum = sum + 0.5*(s1+s2)*dcol(m) s1 = s2 endfor suml = suml + sqrt(sum/uu) endfor sumlor[n] = suml endfor ; -------------------------------------------------------------- ; Caclulate Lorentz opacity ; -------------------------------------------------------------- fac = 2./!pi u = fac*(sumstr/sumlor)^2. fac = 1./(4.*delnu) y = fac*(sumlor^2.)/sumstr f = dblarr(ntau) for n = 0, ntau-1 do begin if(u[n] le 4.d) then begin f[n] = u[n]/sqrt(1.+0.515*u[n]) endif else begin f[n]=(0.5d + 4.d * u[n] - 0.03125/u[n])/sqrt(twopi*u[n])-1.d endelse endfor alor = twopi*y*f ; -------------------------------------------------------------- ; Cacluale Doppler opacity ; -------------------------------------------------------------- fac = 4./sqrt(!pi) w = fac*(sumstr/sumsqr)^2./wdop fac = 1./(4. * delnu) v = fac*wdop*sumsqr^2./sumstr g = dblarr(ntau) for n = 0, ntau-1 do begin if (w[n] le 1.d-3) then begin g[n] = w[n] endif else if (w[n] le 4.) then begin g[n]=3.8035d0*(alog( (w[n]+3.4097d0)/(w[n]+16.5903d0)) + 1.5822d0) endif else begin dum=alog(sqrt(!pi)*w[n]) g(n)=0.1462d0+0.66666667d0*dum^1.5d0+0.3703d0/dum endelse endfor adop = 2.*v*g a = sqrt(adop^2. + alor^2.) if idebug then begin print, 'tmp' print, tmp[0:2], form='(3e11.3)' print, '..', tmp[ntau-2:ntau-1], form='(A11,2e11.3)' print, 'sumstr' print, sumstr[0:2], form='(3e11.3)' print, '..', sumstr[ntau-2:ntau-1], form='(A11,2e11.3)' print, 'sumsqr' print, sumsqr[0:2], form='(3e11.3)' print, '..', sumsqr[ntau-2:ntau-1], form='(A11,2e11.3)' print, 'wdop' print, wdop[0:2], form='(3e11.3)' print, '..', wdop[ntau-2:ntau-1], form='(A11,2e11.3)' print, 'w' print, w[0:2], form='(3e11.3)' print, '..', w[ntau-2:ntau-1], form='(A11,2e11.3)' print, 'g' print, g[0:2], form='(3e11.3)' print, '..', g[ntau-2:ntau-1], form='(A11,2e11.3)' print, 'v' print, v[0:2], form='(3e11.3)' print, '..', v[ntau-2:ntau-1], form='(A11,2e11.3)' print, 'adop' print, adop[0:2], form='(3e11.3)' print, '..', adop[ntau-2:ntau-1], form='(A11,2e11.3)' print, 'alor' print, alor[0:2], form='(3e11.3)' print, '..', alor[ntau-2:ntau-1], form='(A11,2e11.3)' print, 'a' print, a[0:2], form='(3e11.3)' print, '..', a[ntau-2:ntau-1], form='(A11,2e11.3)' end return, a end