NB. Solution of CO rotational lines in optically thick NB. conditions by Coupled Escape Probability method in NB. spherical geometry. NB. Usage: 'R fsb xrs N tau SS Ec' =. RunCO fR; CC; T; N_H NB. fR are the relative shell radii on an arbitrary scale. NB. CC = CO column density in radial direction, T = temperature NB. N_H = proton density / unit volume NB. xrs (x-reshaped): $xrs ->(nshells,nlevels) RunSPCO=: monad define 'fR CC T N_H'=: y mu_out=: int01 50 NB. equal mu steps (51 values) mu_in=: |.(}: cos 0.5p1*int01 50),0 NB. equal angle steps print 'N_H =';N_H;' CO column =';CC;' Temperature =';T SetCO '' NB. read in data for CO molecule xx=. ABF0 '' NB. solve populations w/ preset NRB icount=: 0 NB. set counter to track number of iterations xx=. Iter^:_ xx NB. iterate populations to convergence print 'Single shell escape probability solution: ';icount;' iterations.' print (nlv reshape xx) NB. print level populations of single shell approx. Shells fR NB. set up for shells defined by boundaries "fR" xx=. , ns#"0 xx NB. replicate 1 shell solution for all shells icount=: 0 NB. reset counter xx=. Iter^:_ xx NB. iterate populations to convergence if. 0< +/-.(\:~xx)=xx do. print' * * Population inversion !! * *' end. print 'CEP solution: ';icount;' iterations.' xrs=. nlv reshape xx NB. fractional populations per sublevel N=. N_H*CO_rat*g*"1 xrs NB. actual populations of levels xlmxu=. (-/"1~xrs)+"2 II NB. table of (xl_j - xu_j) xrr=. inz{"1 ,"2 (1|."1 xrs)%xlmxu NB. xrr = xu_j/(xl_j-xu_j) SS=. (3.9729e_16%lam^3)*"1 xrr NB. line source functions kap=. C_tau*"2 (-/"1~xrs) NB. kappas from array of (x_l-x_u) kap=. pdiag"2 }."1 }:"2 kap NB. reformat kappas jemiss=. kap*SS NB. j(emission)= source*kappa Pesc=: 1- +/ |: MMspy NB. prob emission from shell escapes nebula DelVol=. DEL Vol=. (4p1%3)*R^3 NB. physical volume of each shell JemV=. 4p1* DelVol* jemiss NB. line emission in each shell Flx=. +/Pesc* JemV NB. total escaping flux in each line dt=. C_tau*"2 clmn*xlmxu NB. delta taus from array of (x_l-x_u) tau0=. inz3{ ,(0,+/\dt) NB. get line optical depths print 'J''->J |Wavelength|Line Tau |Line Flux' print |:<"0(>uidx;lidx;lam;({:tau0);Flx) R;fR;xrs;N;tau0;SS;Flx ) load '/your_path/splines.ijs' NB. load spline interpolator SetCO=: monad define NB. read and setup data for CO rotational lines CO_rat=: 5e_4 NB. CO/H ratio by number DR=: ({:fR)-({.fR) NB. total fractional shell thicknes ns=: # dfR=: DEL fR NB. ns = # of shells. dfR = fractional shell thicknesses PR=: CC%(CO_rat*N_H) NB. PR = physical radial thickness of whole shell R=: (PR%DR)*(({.fR),({:fR)) NB. shell radii in physical units clmn=: CC NB. clmn = CO column through each shell Elev=: readB'/your_path/CO_rot_E12.dat' NB. energy of levels II=: unit nlv=: #Elev NB. nlv= numb of levels; II= nlv x nlv unit matrix AA=: readB '/your_path/CO_rot_AA12.dat' NB. Einstein A matrix nln=: #inz=: I. 0~:,AA NB. inz = indicies of non-zero elements of AA matrix EE=: |-/~Elev NB. nlv=number of levels; EE=energy difference matrix Eln=: inz{,EE NB. Eln=energy of spectral lines lam2=: 1.239842e_4%EE NB. hypothetical line wavelengths in cm lam=: inz{,lam2 NB. actual lines (transitions with non-zero A's) ordr=: \: lam NB. order by decreasing wavelenth of transition lam=: ordr {lam NB. reordered wavelength list inz=: ordr {inz NB. reorder indices of transitions with non-zero A's ns=: 1 NB. 1st approximation is a single shell inz2=: inz+/~ ,0 inz3=: inz+/~ ,(*:nlv)*(0 1) uidx=: inz{ ,(nlv,nlv)$i. nlv NB. indices of upper levels of transitions lidx=: inz{ ,|:(nlv,nlv)$i. nlv NB. indices of lower levels g=: >: +: i. nlv NB. statistical weights (2*J+1) gc=. |:(# #"_1 ])g NB. nlv x nlv w/ each column a g-value gg=: |:ggt=: (-.lt#g)+(lt#g)*g%/g NB. g_up/g_low (& transpose) for up>low Dnu=: (3211.3* %:T)%lam2 NB. Doppler widths of lines Dnu0=: inz{,Dnu C_tau=: (lam2^2)*AA*gc%8p1*Dnu NB. coefficient for tau (for n_i=N_i/g_i) Exp=: ^_11604.506*EE%T NB. e^(-EE/kT) CC12=. N_H* readB '/your_path/C12_all.dat' NB. CO-H2 collisional rates from Flower and Launay 1985 MNRAS TT=. 10 20 40 60 100 250 NB. TT = temperatures of Flower-Launay tables Cij=: (TT;CC12) splint T NB. interpolate collision de-excitation C_ij=: Cij+(lt nlv)*Exp*(|:Cij) NB. compute collisional excitation rates ) Shells=: monad define NB. setup indicies and columns for multi-shell case R=: (PR%DR)*fR NB. shell radii in physical units ns=: # dfl=. DEL y NB. ns= # of shells; dfl= fractional thickness of shells clmn=: dfl*CC NB. clmn= CO column through each shell print 'Number of shells =';ns;' Total thickness = ';PR;'cm.' inz2=: inz+/~ ,(*:nlv)*i. ns NB. indicies for non-zero elements of higher inz3=: inz+/~ ,(*:nlv)*i. >:ns NB. order (3D) arrays ) NB. reshape=: dyad : '|: y$~x,(x%~#y)' reshape=: [: |: ] $~ [ ,[ %~ [: #] NB. vector to 2-D ABF=: monad define NB. ABF evaluates level population equations yy=. nlv reshape y NB. (ABF outputs error - correct y returns 0) kap=. C_tau*"2 (-/"1~yy) NB. kappas from array of (x_l-x_u) dt=. C_tau*"2 clmn*(-/"1~yy) NB. delta taus from array of (x_l-x_u) dt0=. inz2{ ,dt kap0=. inz2{ ,kap tau0=. inz3{ ,tau=: 0,+/\dt NB. get coupled escape probability matricies MM=. R& SHELL"1 |: kap0 NB. get coupled escape probability matricies MMspy=: MM NB. keep a copy for finding escaping flux pp=. MM p_i yy NB. get all net radiative brackets GM=. C_ij +"2 pp*"2 AA NB. C_ij + p_ij*A_ij 's OUT=. II*"1/~ +/"2 ggt*"2 GM NB. diagonal terms: rates out of levels GM=. g 0}"2 (gg*"2 GM)- OUT NB. 1st row is: Sigma x_i*g_i = 1 B=. 1,(<:nlv)#0 NB. r.h.s. of system ,|: (GM mm"1 yy)-"1 B NB. evaluate equations ) ABF0=: monad define NB. initial populations, no transfer pp0=. 0.5 NB. assume NRB of 1/2 (whatever ...) GM=. (pp0*AA)+C_ij NB. collisions+radiative rates into k OUT=. II* +/ggt*GM NB. rates out of each level k GM=. g 0}(gg*GM)- OUT NB. in(gg*GM)-OUT with 1st sum g*x=1 GM %.~ 1,(<:nlv)#0 NB. solve for populations (linear system) ) Iter=: monad define NB. Newton's method for n equations del=. y*var=. 0.01 NB. fractional variation for num. derivative ax=. ((n,n)$1) +(var*unit n=. #y) dx=. ax*"1 y NB. array of varied x-vectors F0=. ABF y NB. F(x0) DF=. |: (ABF"1 dx)-"1 F0 NB. delta F's Jac=. DF %"1 del NB. (delta F)/delta = Jacobian matrix delta=. (-F0) %. Jac NB. Newton corrections to xx icount=: >: icount NB. increment counter |y+ (y limit delta) NB. apply corrections (no neg. populations) ) limit=: dyad : '(*y)*(0.5*x)<. |y' NB. limit deltas to 50% of y p_i=: dyad define MMij=. x xlmxu=. (-/"1~y)+"2 II NB. add unit diagonal to avoid NaNs xz=. inz{"1 ,"2 (1|."1 y)%xlmxu xij=. |:"2 (|:xz)%/"1 (|:xz) xz=. |: +/"1 xij* MMij NB. Equation (35) frame=. (ns,nlv,nlv)$0 ($frame)$ (,xz)(,inz2)},frame ) NB. Usage R SHELL kap SHELL=: dyad define n=. _1 Rx=. x kap=. y yo=. Rx;kap;mu_out yi=. Rx;kap;mu_in max=. #dr=: DEL Rx Pmat=. (1,max)$0 NB. Initialize (stripped off later) while. max > n=. >:n do. roots=. (n{Rx)+ (((n+1){Rx)-(n{Rx))* Gaussr pp=. (yi&IN"0 roots)+(yo&OT"0 roots) 'fR1 fR2'=. (n,(n+1)){fR dfRn=. fR2-fR1 froot=. fR1+ dfRn* Gaussr NB. roots in fR corresponding to R pp=. pp* *:froot NB. m_ij(r) r^2 Pmat=. Pmat, ((fR2^3)-(fR1^3))%~3* +/pp* dfRn*Gaussw end. NB. multiply by ratio of shell volumes (Vol_j/Vol_i) : R_Vol=. dV %/ dV=. DEL Rx^3 R_kap=. kap %/ kap |:(}. Pmat)* R_Vol*R_kap ) NB. get probability of absorption in each shell for inward rays NB. kap is opacity: d_tau= kap* d_rho NB. Usage: (R;kap;mu) IN r0 --> Probability matrix (R x R) IN=: dyad define 'R kap mu'=. x sth=. %: 1- *:mu NB. sin(theta) for mu=cos(theta) rm2=. *: rmin=. sth* r0=. y NB. rmin is closest to center of mu ray NB. get mumax = mu of rays grazing R's < r0 (0's beyond r0) mumax=. %: 0>. 1- *:R%r0 NB. id=0 if ray with that mu never cuts or touches given R: id=. mumaxr0) 0}"(0 1) id end. z=. %: 0>. (*:R)-/rm2 NB. dist between rmin & R along rays ns=. # Dz=. DEL z NB. ns = # of shells z0=. mu*r0 NB. dist from r0 to rmin point jo=. {.I. r0 Prob(R) OT=: dyad define 'R kap mu'=. x rm2=: y*y*(1- *:mu) NB. rmin^2 (rmin= ray's min dist to center) z=: %: 0>. (*:R)-/rm2 NB. dist between rmin & R along rays z0=: mu*y NB. distance from z to rmin rho=: 0 >. z-"1 z0 NB. distances from z to shell intersections drho=: DEL rho NB. segments of ray cut off by shells dtau=. kap* drho NB. optical depth of shell cuts tau=. 0, +/\ dtau NB. taus from z to shell edges along rays PAB=. -DEL Etas tau NB. fraction reaching segment & absorbed NB. sum over all (+) mu's (1/2 integral from 0 to 1): -: +/"1 (DEL mu)*"1 BIN"1 PAB ) NB. Roots & weights for 10-point Gaussian quadrature over shell thickness: Gaussr=: 0.013046735741414 0.067468316655508 0.16029521585049 0.28330230293538 0.42556283050918 Gaussr=: Gaussr,(1- |.Gaussr) Gaussw=: 0.033335672154344 0.074725674575292 0.10954318125799 0.134633359655 0.14776211235738 Gaussw=: Gaussw, |.Gaussw 'logt leta Dleta'=: |: mread'/your_path/Log_eta_table.dat' Etas=: monad define pick=. I. 1>t=. ,y NB. arguments less than unity eta=. 10^ (logt;leta;Dleta) spltrp 10^. t>. 0.1 small=. ps_eta pick{t NB. power series for tau<1 ($y)$ small pick} eta NB. insert "small" & restore dimensions ) NB. power-series approximation to Eta for small tau. NB. Usage: "ps_coef n" gives 1st n coefficients of power series ps_coef=: monad : '(_1p_0.5^n)%(!n)*(%:1+n=. i.y)' ps_eta=: (ps_coef 14)&p. f. NB. Usage: Eta=. ps_eta tau NB. 13-place accuracy for tau<1. pdiag=: (<0 1)&|: NB. get principle diagonal of matrix unit=: =@i. NB. unit N ==> N x N unit matrix lt=: >:/~@i. NB. generate lower triangular matrix ref_axis=: monad : '(}: -|. y), y' NB. reflect x-axis about 0 ref_data=: monad : '(}:"1 |."1 y),"1 y' NB. reflect data points about 0 NB. tau1=: 0.01 NB. first optical depth point NB. npd=: 5 NB. number of points/decade NB. T=: 4500 NB. Half-depth of slab tau_grid=: monad define 'tau1 npd T'=: y T,~(i. +/tJ |Wavelength|Line Tau |Line Flux NB. ┌──┬──┬─────────┬────────┬──────────┐ NB. │1 │0 │0.260076 │7.87508 │1.68372e20│ NB. ├──┼──┼─────────┼────────┼──────────┤ NB. │2 │1 │0.130038 │26.6188 │1.82802e21│ NB. ├──┼──┼─────────┼────────┼──────────┤ NB. │3 │2 │0.0866919│47.7926 │5.55719e21│ NB. ├──┼──┼─────────┼────────┼──────────┤ NB. │4 │3 │0.0650191│60.7554 │9.85644e21│ NB. ├──┼──┼─────────┼────────┼──────────┤ NB. │5 │4 │0.0520153│57.5402 │1.29697e22│ NB. ├──┼──┼─────────┼────────┼──────────┤ NB. │6 │5 │0.043346 │40.7839 │1.41164e22│ NB. ├──┼──┼─────────┼────────┼──────────┤ NB. │7 │6 │0.0371537│23.8456 │1.25673e22│ NB. ├──┼──┼─────────┼────────┼──────────┤ NB. │8 │7 │0.0325095│11.8669 │9.31927e21│ NB. ├──┼──┼─────────┼────────┼──────────┤ NB. │9 │8 │0.0288973│5.49239 │5.64365e21│ NB. ├──┼──┼─────────┼────────┼──────────┤ NB. │10│9 │0.0260076│2.2607 │2.9044e21 │ NB. ├──┼──┼─────────┼────────┼──────────┤ NB. │11│10│0.0236433│0.850131│1.27715e21│ NB. └──┴──┴─────────┴────────┴──────────┘