load '/your_path/read_beams45.ijs' DTR=: 1p1%180 NB. convert degrees to radians NB. Compute emergent Stokes parameters for array of longitudes NB. latitudes over illuminated face of planet for given phase. RUN_PLANET=: monad define 'albedo phase'=: y 'muin mubin abin Su'=: Beam_read45 albedo NB. Su are "Stokes parameters" from running "Pol_beam.ijs". BUT they are NB. not real intensities: they are photon fluxes per unit surface area, NB. *not* per unit *projected* area as seen from angle "mu". Thus we must NB. divide the Su's by mu(=cos theta) convert them to specific intensities. Su=: mubin%"2~ Su NB. convert photon numbers/area to specific intensity Su=: (#mubin)*(1p1%~#abin)*Su NB. also, remove included (del mu)*(del phi) NB. ------------------ test of Lambert surface --------------------------- NB. n=. $(|:"_1 Su) NB. rearranged dimension of Su NB. Set intensity of scattered radiation (I,U,Q)=(1,0,0) for all angles: NB. Su=: |:"_1 (n$1 0 0) NB. and restore dimension of original NB. ---------------------------------------------------------------------- NB. 'fret' is the grid of mu control points for the least-squares fits: fret=. 0 0.05 0.15 0.3 0.5 0.7 0.85 0.95 1.0 ALLOUT=. (fret;muin) ls_splines (|: ,"3 Su) NB. construct l.sq. spline fits low=. phase - 90 NB. longitude integration runs from (-90+phase)-->(+90) high=. 90 del=. high - low ll=: DTR* }.(low- -:del%nl)+ del*int01 nl=. 41 pp=: DTR* }.(_45%np)+ 90*int01 np=. 29 phase=: DTR* phase wxz=. weights '' NB. projected areas of longitude-latitude cells EE=. (np,nl)$0 NB. placeholder for mu_in(phi,lon) table n1=. _1 while. (n1=. >:n1)< nl do. l=. n1{ll n2=. _1 while. (n2=. >:n2)< np do. ph=. n2{pp E=. (cos l)* cos ph NB. mu_in for illuminating rays EE=. E (:n1)< nl do. l=. n1{ll n2=. _1 while. (n2=. >:n2)< np do. ph=. n2{pp NB. n1{n2{ SX is the interpolated table for this mu_in(l,ph) S=. (n1{n2{ SX) PLANET phase, l, ph SS=. S ( (I,Q,U) PLANET=: dyad define 'phase lon phi'=. y cphi=. cos phi th_i=. acos cth_i=. cphi* cos lon sth_i=. %: 1- *: cth_i th_o=. acos cth_o=. cphi* cos(phase-lon) sth_o=. %: 1- *: cth_o cph_o=. ((cos phase)-cth_i*cth_o)%(sth_i*sth_o) ph_o=. acos cph_o 'I Q U'=. x Scat cth_i, th_o, ph_o alph=. acos (sin phi)%sth_o s2a=. %: 1- *:c2a=. cos +:alph if. lon>phase do. s2a=. -s2a end. NB. rotate Q,U to meridian plane: I,((Q*c2a)+(U*s2a)),((-Q*s2a)+(U*c2a)) ) NB. Linear interpolation in S tables for the NB. specific theta(out) and phi(out): Scat=: dyad define 'mu_in th_o ph_o'=. y mu_o=. cos th_o n=. mubin IdX mu_o SS1=. n{"2 x SS2=. (>:n){"2 x 'mu1 mu2'=. (n,>:n){ mubin del=. mu2- mu1 g2=. del%~ mu_o-mu1 g1=. del%~ mu2-mu_o SS=. (g1*SS1)+ g2*SS2 n=. abin IdX ph_o SS1=. n{"1 SS SS2=. (>:n){"1 SS 'ph1 ph2'=. (n,>:n){ abin del=. ph2 - ph1 g2=. del%~ ph_o-ph1 g1=. del%~ ph2-ph_o (g1*SS1)+ g2*SS2 ) weights=: monad define pv=: 0,({. pp)+pp NB. boundaries in latitude h=. -: {. DEL ll NB. 1/2 the longitude step lv=: (ll-h), h+{:ll NB. boundaries in longitude NB. area element dA = (cos lat)* d(lat)*d(long) NB. projection factor for dA= (cos lat)*cos(long-phase) pq=. (*: cos pp)* DEL pv lq=. (cos ll-phase)* DEL lv pq*/lq NB. construct 2-D table ) NB. x IdX x0 --> index where x0 fits in x-array, e.g.: NB. 0 1 2 3 4 5 6 <- x values NB. | | | | | | | NB. <-- 0 | 1 | 2 | 3 | 4 | 5 ---> IdX values IdX=: ] I.~ [:}. [:}: [ NB. (}.}: t) I. x0 NB. Routine to find cubic spline fit to data by least squares. NB. t= boundaries of cubic segments, (xi,yi) = data points. NB. Usage: 't z zp'=. OUT=. (t,xi) ls_spline yi NB. z and zp are values of spline fit f(x) and df/dx at t's. NB. (z & zp are sufficient to evaluate the fit f(x) at any x.) NB. This version assumes multiple "yi" data arrays ls_splines=: dyad define 'tt xi'=. x xi=. (o=. /:xi){xi yi=. o{"1 y NB. sort xi into an ascending sequence idx=. tt IdX xi NB. idx is t interval in which each xi resides NB. if first or last interval(s) are empty of data, delete them t=. ({. idx)}.((2-#tt)+{:idx)}. tt (t;xi) ls_funct"1 yi NB. do each yi one at a time ) NB. Use least squares to find parameters of spline fits. NB. (see description on webpage under "Maths") ls_funct=: dyad define 't xi'=. x yi=. y Np1=. 2+ Nm1=. <: N=. #h=. (}.-}:)t bx=. <"1 |: >(i. n=. #xi); t IdX xi v=. _1+ u=. bx{ h%"1~ xi -/ }:t uu=. *: u=. u bx}0$~n,N vv=. *: v=. v bx}0$~n,N aa=. +/ *:al=. (1+2*u)*vv ab=. +/ al*bt=. uu*(1-2*v) ag=. +/ al*gm=. h*"1 u*vv ad=. +/ al*dl=. h*"1 uu*v bd=. +/ bt*dl [bg=. +/ bt*gm [bb=. +/ *:bt dd=. +/ *:dl [gd=. +/ gm*dl [gg=. +/ *:gm D=. +:(0,~ +/ yi*al)+(0, +/ yi*bt) G=. +: (0,~ +/ yi*gm)+(0, +/ yi*dl) bx11=. >: &.> bx00=. ;/ ,.~i. N bx01=. (0 1)&+ &.> bx00 bx10=. (1 0)&+ &.> bx00 A=. (bb bx11}E) + aa bx00}E=. 0$~,~Np1 A=. +: (ab,ab) (bx01,bx10)}A B=. (ag bx00}E) + bd bx11}E B=. +: (bg,ad) (bx01,bx10)}B E=. (gg bx00}E) + dd bx11}E E=. +: (gd,gd) (bx01,bx10)}E bx00=. }: bx00 [ bx01=. }: bx01 bx02=. (0 2)&+ &.> bx00 hip=. }.h [hi=. }:h C=. (3*hip%hi) bx00}F=. 0$~Nm1,Np1 C=. (3*(hi%hip)-hip%hi) bx01}C C=. (_3*hi%hip) bx02}C F=. hip bx00}F F=. (+: hi+hip) bx01}F F=. hi bx02}F NB. set up matrix of linear system to be solved: M=. (A,.(|:B),.|:C),(B,.E,.|:F),C,.F RHS=. (D,G,Nm1#0) NB. right hand side of system ZZ=. RHS %. M NB. solve for f(x_i) & [df/dx]_i z=. (i. Np1){ZZ >t; z; zp=. (Np1+i. Np1){ZZ ) NB. Least squares spline evaluation using output "OUT" NB. from ls_splines with multiple "yi" data sets. NB. Usage: y eval_spline"2 OUT --> f(y) NB. The argument y may be an array of x-values. eval_splines=: dyad define 't z zp'=. y NB. Array t holds the boundaries of the spline pieces, z is NB. the array of the values of the cubic splines at the NB. t points, and zp (z') the derivatives at those points. N=. # h=. (}.-}:) t v=. _1+ u=. h%~ |: x -/ }:t a=. (1+2*u)*v*v* }:z b=. u*u*(1-2*v)* }.z c=. h*u*v*v* }:zp d=. h*u*u*v* }.zp box=. <"1 |: >(i. #x); t IdX x box { ((#x),N)$,|: a+b+c+d ) NB. Fractional polarization from (I,Q,U) vector or (3,n) array NB. added 1<. to POL=: {. %"1~ [: L }. so result can't exceed unity. POL=: 1<. {. %"1~ [: L }. L=: +/ &. (*:"_) NB. Gives the lengths of vector arrays NB. Angle of plane of polarization wrt reference plane of NB. (I,Q,U) vector or (3,n) array CHI=: [: -: [: atan2r 2 1 { ] atan2r=: 13 : '2p1|((0>1{y)*1p1)+ 2p1+ atan %/y=. Re y' NB. I, Fractional polarization & angle from (I,Q,U) vector (or array) FracPol=: monad define pol=. POL y chi=. CHI y (0{y), pol,: (180%1p1)*chi )