Function Pn,z,zeta
;note that zeta is always set to zero for all the calculations i'm doing, it's just easier to specify that in another program.
pi = 3.14159265
h = 100000
n1 = randomn(seed,h)
n2 = randomn(seed,h)
n3 = randomn(seed,h)
n4 = randomn(seed,h)
n5 = randomn(seed,h)
restore,filename = 'one_mg.sav', /ver
psd = psd_iso(mg, /arc,nres = 2000)
;since the mdi database has gone down, Lucas has given me a quiet sun magnetogram to work with
l = 10. ; this variable makes it easier to edit the elements of c in the main equation
z = findgen(500)/l ; height above photosphere
n = n_elements(psd.k)
binsize = psd.k[2]-psd.k[1]
A = 0
for i = 0,n-1 do begin $
A = A + binsize*psd.k[i]*psd.psd[i]*exp(-2.*psd.k[i]*z) & $
endfor
A = A*2.*pi
B = 0
for i = 0,n-1 do begin $
B = B + binsize*(psd.k[i])^3.*psd.psd[i]*exp(-2.*psd.k[i]*z) & $
endfor
B = B*2.*pi
D = 0
for i = 0,n-1 do begin $
D = D + binsize*(psd.k[i])^2.*psd.psd[i]*exp(-2.*psd.k[i]*z) & $
endfor
D = D*2.*pi
q = (b/a)^0.5
sigma = (a)^0.5
c = D/((a*b))^0.5
; these are the same spectral integrals as in my program integrateS, because the function doesn't remember them
G = exp(-0.5*(zeta)^2)/(4*sqrt(2)*(pi)^(1.5))*(((sqrt(1-c[z*l]^2)*n5 + c[z*l]*zeta)^3) + ( (1-c[z*l]^2)*((n3^2 + n4^2) -0.5*(n1^2 +n2^2)*(sqrt(1-c[z*l]^2)*n5 + c[z*l]*zeta)) ) + (1-c[z*l]^2)*(sqrt(2)*n2*n3*n4 + (1/sqrt(2))*(n3^2 - n4^2)*n1)) * q[z*l]^3
G = abs(G)
return, mean(G)
; this is the full expression for null densities