NB. Solution of [O I] fine structure lines in optically thick NB. conditions by Coupled Escape Probability method. NB. To start a new (N) calculation: NB. output=. 'N' Run fsb; 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' Run=: dyad define Setup y if. (x='B')+. (x='N') do. NB. if New or Back if. x='B' do. NB. if Back read in stored xx xx=. readB '/your_path/xx_back.dat' elseif. x='N' do. NB. if New xx=. ABF0 '' NB. solve populations w/ preset p_ij end. else. NB. otherwise approximate x2 & x3 provided xx=. ,> x NB. concatinate x2,x3 end. xx writeB '/your_path/xx_back.dat' xx=. Iter^:_ xx NB. iterate populations to convergence 'x2 x3'=. box2 xx NB. divide solution vector into x2 & x3 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)) E21=. Dnu21* +/(Esc tau21)*S21 NB. line cooling = (intensity/Dnu)*Dnu E32=. Dnu32* +/(Esc tau32)*S32 NB. final tau's set in last call to ABF print 'total tau(2->1)=';(ns{tau21);'; total tau(1->0)=';(ns{tau32) string=. 'Line emission: 63.2 mic=' print string;E21;' 145.5 mic=';E32;' Ratio =';(E32%E21) print 'Emission /O atom/cc: (63m)';(E21%CC);' (145m)';(E32%CC) x1;x2;x3;N1;N2;N3;S21;S32;E21;E32 ) Setup=: monad define NB. setup data for [O I] fine structure lines 'fsb CC T N_H'=: y NB. fsb = slab boundaries on [0,1] interval OH_rat=: 5e_4 NB. O/H ratio by number print 'N_H =';N_H;' Oxygen column =';CC;' Temperature =';T ns=: # dfl=. DEL fsb NB. ns = # of slabs. CC= total column of oxygen LL=: CC%(OH_rat*N_H) NB. LL = total slab thickness clmn=: DEL fsb*CC 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 slabs =';ns;' Total thickness = ';LL;'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 NB. evaluate errors in current x2 & x3 values 'x2 x3'=. box2 y NB. split input vector into x2 & x3 parts x1=. g1%~ 1-(g2*x2)+(g3*x3) M21=. MMs tau21=: 0, +/\dt21=: O_21*clmn*(x1-x2) NB. line optical depths M32=. MMs tau32=: 0, +/\dt32=: O_32*clmn*(x2-x3) NB. and escape matrices p21=. (dt21;M21) p_i x2;x1 NB. compute net radiative brackets p32=. (dt32;M32) p_i x3;x2 NB. the rate equations are a11*x2 + a12*x3 - rhs1 = 0 NB. and a21*x2 + a22*x3 - rhs2 = 0 NB. Note: We have used "g1*x1+g2*x2+g3*x3=1" to explicitly eliminate x1. a11=. -(A21*p21)+C21+(C32*Exp32%g2)+(g2%g1)*C21*Exp21 a12=. (A32*p32%g2)+(C32%g2)-C21*Exp21%g1 a21=. ns# (-C32*Exp32)+(g2%g1)*C31*Exp31 a22=. (A32*p32)+C32+C31*(1+Exp31%g1) rhs1=. ns# -C21*Exp21%g1 rhs2=. ns# C31*Exp31%g1 A1=. (x2*a11)+(x3*a12) A2=. (x2*a21)+(x3*a22) (A1-rhs1),(A2-rhs2) NB. these are the errors in x2 & x3 ) ABF0=: monad define NB. initial populations (no radiative transfer) 'p21 p32'=. 0.3 0.9 NB. fraction assumed escaping a11=. -((A21*p21)+C21+(C32*Exp32%g2)+(g2%g1)*C21*Exp21) a12=. ((A32*p32%g2)+(C32%g2)+(-C21*Exp21%g1)) a21=. ((-C32*Exp32)+(g2%g1)*C31*Exp31) a22=. (C31+(A32*p32)+C32+C31*Exp31%g1) rhs=. g1%~ (-C21*Exp21),(C31*Exp31) 'xx2 xx3'=. rhs %. >(a11,a12);(a21,a22) 'x2 x3'=. box2 xx=. (ns&#) xx2,xx3 NB. replicate for all slabs x1=. g1%~ 1-(g2*x2)+(g3*x3) tau21=: 0, +/\dt21=: O_21*clmn*(x1-x2) NB. evaluate line optical depths tau32=: 0, +/\dt32=: O_32*clmn*(x2-x3) xx ) Iter=: monad define NB. Newton's method for n equations del=. y* var=. 0.0002 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 |y+ (y limit delta) NB. apply corrections (no neg. populations) ) limit=: dyad : '(*y)*(0.1*x)<. |y' NB. limit deltas to 10% box2=: [: <"1 ] $~ 2,2 %~ # NB. split concatinated vector NB. Calculate the net radiative brackets, i.e., the NB. p^i values of Elitzur & Ramos (2006), eqn (35) NB. Usage: (tau vector; Mij matrix) p_i (x_u;x_l) p_i=: dyad define 'dtau 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) dtau%~ +/"1 xij*MMij NB. Equation (35) ) NB. Calculate the M_ij matrix, Equation (18) NB. Elitzur & Ramos (2006) MNRAS 365, 779. NB. Usage: MMs tau (tau array, not dtau). MMs=: monad define z=. - -: Alphas |-/~y (}.}."1 z)+(-}:}."1 z)+(-}.}:"1 z)+ }:}:"1 z ) load'/your_path/splines.ijs' 'logt lalpha Dlal'=: readB '/your_path/log_alpha_table_161.dat' Alphas=: monad : '10^ (logt;lalpha;Dlal) spltrp 10^. y' unit=: =@i. NB. unit N ==> N x N unit matrix 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=: T,~(i. +/t2)=5.64224e_10 (0->2)=4.80021e_10 (0->1)=8.16784e_10 NB. total tau(2->1)=40.435; total tau(1->0)=7.31935 NB. Line emission: 63.2 mic=0.243511 145.5 mic=0.0237171 Ratio =0.0973965 NB. Emission /O atom/cc: (63m) 2.43511e_20 (145m) 2.37171e_21 Esc=: monad : '2p1* (|. Dal)+ Dal=. DEL Alphas y' NB. Jbar=. +/(Esc tau)*S(tau) NB. Calculate the profiles of the lines: NB. load '/home/jph/J6/CEP_PP/Flux.ijs' NB. load '/home/jph/J6/CEP_PP/Alpha.ijs' NB. xx=. 5*int01 100 NB. p21=. tau21 FluxF S21 NB. p32=. tau32 FluxF S32 NB. plot the profiles: NB. load'plot' NB. 'pensize 3' plot ((}: -|.xx),xx);>((}:|.p21),p21);((}:|.p32),p32)