NB. Solution of [O I] fine structure lines in optically thick NB. spherical shells by the Coupled Escape Probability method. NB. To start a new (N) calculation: NB. output=. 'N' Run R; CC; T; N_H NB. To use previous populations x2 & x3 as 1st guess: NB. output=. (x2;x3) Run fsb; CC; T; N_H NB. If run crashes, to step back (B) to previous good run: NB. output=. 'B' Run fsb; CC; T; N_H NB. output = 'x1 x2 x3 N1 N2 N3 S21 S32 E21 E32' SP_OI=: dyad define Setup y mu_out=: int01 60 NB. equal mu steps (61 values) mu_in=: |.(}: cos 0.5p1*int01 60),0 NB. equal angle steps if. (x='B')+. (x='N') do. NB. if New or Back if. x='B' do. xx=. readB '/your_path/xx_back.dat' elseif. x='N' do. xx=. ABF0 '' NB. solve populations w/ preset p_ij end. else. xx=. ,> x NB. first guess at x2;x3 end. xx writeB '/your_path/xx_back.dat' xx=. Iter^:_ xx NB. iterate populations to convergence 'x2 x3'=. box2 xx print' Iterations converged !' x1=. g1%~ 1-(g2*x2)+(g3*x3) 'N1 N2 N3'=: N_H*OH_rat*(g1,g2,g3)* >x1;x2;x3 S21=. (3.9729e_16%lam21^3)*x2%(x1-x2) NB. line source functions: S32=. (3.9729e_16%lam32^3)*x3%(x2-x3) NB. (2ch/lam^3)*(n_u/(n_l-n_u)) print 'total tau(2->1)=';(ns{tau21);'; total tau(1->0)=';(ns{tau32) x1;x2;x3;N1;N2;N3;S21;S32 ) Setup=: monad define NB. setup data for [O I] fine structure lines 'fR CC T N_H'=: y NB. fR = fractional radii of spherical shells OH_rat=: 5e_4 NB. O/H ratio by number print 'N_H =';N_H;' Oxygen column =';CC;' Temperature =';T ns=: # dfR=: DEL fR NB. ns = # of shells. CC= total column of oxygen DR=: ({:fR)-({.fR) NB. total fractional shell thickness PR=: CC%(OH_rat*N_H) NB. PR = physical radius of sphere R=: (PR%DR)*fR NB. shell radii in physical units clmn=: CC*dfR%DR NB. clmn = oxygen column through each slab 'lam21 lam32'=: 63.184e_4 145.53e_4 NB. line wavelengths in cm 'g1 g2 g3'=: 5 3 1 NB. statistical weights 'A21 A31 A32'=: 8.92e_5, 0, 1.74e_5 NB. A31 actually 1e_10, neglected Dnu21=: (3211.3* %:T)%lam21 NB. Doppler width of line Dnu32=: (3211.3* %:T)%lam32 O_21=: (lam21^2)*A21*g2%8p1*Dnu21 NB. coefficient for tau_21 (based on n_1=N_1/g_1) max_t=. O_21*CC%g1 NB. tau21 if all in J=2 level (CC divided by g_1=5) O_32=: (lam32^2)*A32*g3%8p1*Dnu32 NB. coefficient for tau_32 (based on n_2=N_2/g_2) print 'Number of shells =';ns;' Shell thickness = ';PR;'cm. Max tau21=';max_t ee21=. 1.43878% lam21 NB. energy of transition/k Exp21=: ^-ee21% T NB. e^(-E_21/kT) ee32=. 1.43878% lam32 Exp32=: ^-ee32% T Exp31=: ^-(ee21+ee32)% T NB. data for hydrogen collisions: Abrahamsson et al.(2007) ApJ 654,1171 'Tval OC21 OC31 OC32'=. readB '/your_path/OI_H_fine_coll.dat' C21=: 1e_11*((Tval;OC21) splint T) NB. interpolate [O I] collision rates C31=: 1e_11*((Tval;OC31) splint T) NB. with H, tabulated for 502)=';C21;' (0->2)=';C31;' (0->1)=';C32 'C21 C31 C32'=: N_H* C21, C31, C32 NB. multiply rates by H number density ) ABF=: monad define 'x2 x3'=. box2 y x1=. g1%~ 1-(g2*x2)+(g3*x3) tau21=: 0,+/\ dt21=: clmn* kap21=: O_21*(x1-x2) NB. line opacities tau32=: 0,+/\ dt32=: clmn* kap32=: O_32*(x2-x3) M21=: SHELL R;kap21 NB. coupling matrix for 2-1 line M32=: SHELL R;kap32 NB. and for 3-2 line p21=. M21 p_i x2;x1 NB. compute net radiative brackets p32=. M32 p_i x3;x2 NB. the rate equations are a11*x2 + a12*x3 = rhs1 NB. and a21*x2 + a22*x3 = rhs2 NB. Note: We have used "g1*x1+g2*x2+g3*x3=1" to explicitly eliminate x1. U=. unit ns a11=. -U*(A21*p21)+C21+(C32*Exp32%g2)+(g2%g1)*C21*Exp21 a12=. U*(A32*p32%g2)+(C32%g2)-C21*Exp21%g1 a21=. U*(-C32*Exp32)+(g2%g1)*C31*Exp31 a22=. U* C31+(A32*p32)+C32+C31*Exp31%g1 rhs1=. ns# -C21*Exp21%g1 rhs2=. ns# C31*Exp31%g1 AA=. (a11,"1 a12),"2 a21,"1 a22 NB. grand coefficient matrix B=. rhs1,rhs2 NB. constant terms of system NB. AA*(x2,x3)=B; for approximate x2,x3: AA*(x2,x3)-B = error (AA;B) FF y NB. evaluate equations for populations y ) ABF0=: monad define NB. get populations with no radiative transfer 'p21 p32'=. 0.3 0.9 NB. fraction assumed escaping a11=. -((A21*p21)+C21+(C32*Exp32%g2)+(g2%g1)*C21*Exp21)*unit ns a12=. ((A32*p32%g2)+(C32%g2)+(-C21*Exp21%g1))*unit ns a21=. ((-C32*Exp32)+(g2%g1)*C31*Exp31)*unit ns a22=. (C31+(A32*p32)+C32+C31*Exp31%g1)*unit ns rhs1=. ns# -C21*Exp21%g1 rhs2=. ns# C31*Exp31%g1 'x2 x3'=. box2 xx=. (rhs1, rhs2)%. (a11,"1 a12),"2 a21,"1 a22 x1=. g1%~ 1-(g2*x2)+(g3*x3) kap21=: O_21*(x1-x2) NB. evaluate line optical depths kap32=: O_32*(x2-x3) xx ) Iter=: 3 : 0 NB. Newton's method for n equations del=. y* var=. 0.002 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) vF=. ABF"1 dx NB. F at varied x's DF=. |: vF-"1 F0 NB. delta F's Jac=. DF %"1 del NB. (delta F)/delta = Jacobian matrix delta=: -F0 %. Jac NB. Newton corrections to xx print' deltas: ',":delta NB. so we can watch the convergence |y+ (y limit delta) NB. apply corrections (no neg. populations) ) FF=: 4 : '((>0{x) mm y) - (>1{x)' limit=: 4 : '(*y)*(0.1*x)<. |y' NB. limit deltas to 10% box2=: [: <"1 ] $~ 2,2 %~ # NB. Calculate the net radiative brackets, i.e., the NB. p^i values of Elitzur & Ramos (2006), eqn (35) NB. Usage: (coupling matrix MM) p_i (x_u;x_l) p_i=: dyad define MMij=. x 'xu xl'=. y xx=. xu%(xl-xu) NB. xx = xu_j/(xl_j-xu_j) xij=. |: xx%/xx NB. xu_ij=(xu_j/xu_i)*(xl_i-xu_i)/(xl_j-xu_j) +/"1 xij*MMij NB. p_i ) pdiag=: (<0 1)&|: NB. get principle diagonal of matrix unit=: =@i. NB. unit N ==> N x N unit matrix NB. Compute coupling of shells. Usage: SHELL R;kap SHELL=: monad define n=. _1 'Rx 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)) vmap Gaussr pp=. (yi&IN"0 roots)+(yo&OUT"0 roots) Pmat=. Pmat, +/pp*Gaussw NB. integrate over r_i 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. For quadrature roots over interval [0,1], find radii of NB. the shells that divide the *volume* according to the roots vmap=: dyad define 'Ri Rip'=. x vi=. Ri^3 DelV=. (Rip^3) - vi vv=. vi+y*DelV vv^%3 ) 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) OUT=: 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 ) load'/your_path/splines.ijs' NB. read in table of log(tau), log(eta) & 2nd derivatives: 'logt leta Dleta'=: |: mread'/your_path/Log_eta_table.dat' 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 NB. Interpolate in table to find eta(tau)s: 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. 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=: 3 : 0 'tau1 npd T'=: y T,~(i. +/t.(*:r)-/*:(0,r) NB. r_min to shell cuts zz=. }:"1 DEL z NB. incremental distances rat=. zz% {."1 zz NB. (slant/radial) distances dts=. rat* DEL tau NB. tau increments along rays taus=. 0, +/\ dts NB. tau along rays taur=. ({:taus)-"1 taus NB. reverse taus: we're looking in aa=. ^-phi0* taur NB. attenuated contribution: exp(-tau_x) near=. +/s* DEL aa NB. multiply by source fun. and sum up tf=. ({:taus)+"1 taus NB. taus of rays to shells on far side bb=. ^-phi0*tf NB. contrib. from far side near+ +/-s* DEL bb NB. near + far contributions )