load'/your_path/MTs.ijs' load'/your_path/IQMs.ijs' load'/your_path/IQs.ijs' NB. Usage: (mu out array) U_sp (t1, npd, T) NB. where T is the half-thickness of the slab U_sp=: 4 : 0 't1 npd T'=. y ntau=. # tau=. tau_grid2 y mu=. x NB. source term is uniform throughout slab Ss=. ntau# %8p1*T NB. S* = 1/(4 pi (2*T)) F=. 2 MTs tau F4=. 4 MTs tau LL=. 1 MTs tau MM=. 3 MTs tau NN=. 5 MTs tau Id=. =i. ntau NB. identity matrix B=. lam=. 0 NB. placeholders for future lam & B a=. Id - (1-lam)*LL b=. _3%~(1-lam)*MM c=. (_3%8)*(1-lam)*MM d=. Id - (3%8)*(1-lam)*NN Q=. (a,"1 b),"2 c,"1 d NB. form grand coefficient matrix Q=[a,b/c,d] rhs=. (ntau#0),~ Ss + lam*(ntau# B) 's p'=. (2,ntau)$rhs %. Q NB. solve linear eqns: rhs %. (coeff matrix) Iout=. IQMs tau;mu NB. matrix for emergent intensity, etc. 'Itop Q pol'=. (Iout;mu) IQs s;p Flx=: (F mm s)+(F4 mm p) NB. flux array [(+) =up, (-) =down] Fs=. +: 1p1* {.Flx NB. flux scattered out top + bottom Itop;Q;pol;s;p;Fs ) NB. Construct logarithmic grid over the interval [0,HT], NB. symmetric about central palne at tau=HT. NB. tau1 <- first optical depth point, e.g. tau1=: 0.01 NB. npd <- number of points/decade, e.g. npd=: 5 NB. HT <- Half-depth of slab, e.g. HT=: 2 NB. Usage: tau_grid2 tau1 npd HT --> tau grid tau_grid2=: 3 : 0 'tau1 npd HT'=: y T=. +:HT t=. HT,~(i. +/t