;+
; Contains the magnetic_field procedure
; and associated functions
;
; :Author:
; Baptiste Cecconi, Laurent Lamy, Philippe Zarka
;
; :History:
; 2008/02/20: Created
;
; 2008/02/20: Last Edit
;-
;
;-------------------------------------------------------------------------------
function P10,H
return,COS(H)
end
;-------------------------------------------------------------------------------
function P11,H
return,SIN(H)
end
;-------------------------------------------------------------------------------
function P20,H
return,((3.d0*COS(2.d0*H)+1.d0)/4.d0)
end
;-------------------------------------------------------------------------------
function P21,H
return,(sqrt(3.d0) / 2.d0 * SIN(2.d0 * H))
end
;-------------------------------------------------------------------------------
function P22,H
return,(sqrt(3.d0) / 4.d0 * (1.d0 - COS(2.d0 * H)))
end
;-------------------------------------------------------------------------------
function P30,H
return,(5.d0*COS(3.d0*H)+3.d0*COS(H))/8.d0
end
;-------------------------------------------------------------------------------
function P31,H
return,sqrt(3.d0) / 8.d0 / sqrt(2.d0) * (5.d0 * SIN(3.d0 * H) + SIN(H))
end
;-------------------------------------------------------------------------------
function P32,H
return,(sqrt(15.d0) / 8.d0 * (COS(H) - COS(3.d0 * H)))
end
;-------------------------------------------------------------------------------
function P33,H
return,(sqrt(5.d0) / 8.d0 / sqrt(2.d0) * (3.d0 * SIN(H) - SIN(3.d0 * H)))
end
;-------------------------------------------------------------------------------
function P40,H
return,((35.d0 * COS(4.d0 * H) + 20.d0 * COS(2.d0 * H) + 9.d0) / 64.d0)
end
;-------------------------------------------------------------------------------
function P41,H
return,(sqrt(5.d0) / 16.d0 / sqrt(2.d0) * (7.d0 * SIN(4.d0 * H) + 2.d0 * SIN(2.d0 * H)))
end
;-------------------------------------------------------------------------------
function P42,H
return,(sqrt(5.d0)/32.d0*(3.d0+4.d0*COS(2.d0*H)-7.d0*COS(4.d0*H)))
end
;-------------------------------------------------------------------------------
function P43,H
return,(sqrt(35.d0) / 16.d0 / sqrt(2.d0) * (2.d0 * SIN(2 * H) - SIN (4 * H)))
end
;-------------------------------------------------------------------------------
function P44,H
return,(sqrt(35.d0) / 64.d0 * (COS(4.d0 * H) - 4.d0 * COS(2.d0 * H) + 3.d0))
end
;;-------------------------------------------------------------------------------
;function P50,H
;return,(1.d0/128.d0*(63.d0*COS(5.d0*H)+35.d0*COS(3.d0*H)+30.d0*COS(H)))
;end
;;-------------------------------------------------------------------------------
;function P51,H
;return,(sqrt(15.d0)/128.d0*(21.d0*SIN(5.d0*H)+7.d0*SIN(3.d0*H)+2.d0*SIN(H)))
;end
;;-------------------------------------------------------------------------------
;function P52,H
;return,(sqrt(3.d0)*sqrt(35.d0)/64.d0*(-3.d0*COS(5.d0*H)+COS(3.d0*H)+2.d0*COS(H)))
;end
;;-------------------------------------------------------------------------------
;function P53,H
;return,(sqrt(35.d0)/128.d0/sqrt(2.d0)*(-9.d0*SIN(5.d0*H)+13.d0*SIN(3.d0*H)+6.d0*SIN(H)))
;end
;;-------------------------------------------------------------------------------
;function P54,H
;return,(3.d0*sqrt(35.d0)/128.d0*(COS(5.d0*H)-3.d0*COS(3.d0*H)+2.d0*COS(H)))
;end
;;-------------------------------------------------------------------------------
;function P55,H
;return,(3.d0*sqrt(7.d0)/128.d0/sqrt(2.d0)*(SIN(5.d0*H)-5.d0*SIN(3.d0*H)+10.d0*SIN(H)))
;end
;;-------------------------------------------------------------------------------
function DP10,H
return, -sin(H)
end
;-------------------------------------------------------------------------------
function DP11,H
return,COS(H)
end
;-------------------------------------------------------------------------------
function DP20,H
return, (-3.d0/2.d0*sin(2.d0*H))
end
;-------------------------------------------------------------------------------
function DP21,H
return,(sqrt(3.d0) * COS(2.d0 * H))
end
;-------------------------------------------------------------------------------
function DP22,H
return,(sqrt(3.d0) / 2.d0 * SIN(2.d0 * H))
end
;-------------------------------------------------------------------------------
function DP30,H
return, (-3.d0/8.d0*(5.d0*sin(3.d0*H)+sin(H)))
end
;-------------------------------------------------------------------------------
function DP31,H
return,(sqrt(3.d0) / 8.d0 / sqrt(2.d0 ) * (15.d0 * COS(3.d0 * H) + COS(H)))
end
;-------------------------------------------------------------------------------
function DP32,H
return,(sqrt(15.d0) / 8.d0 * (3.d0 * SIN(3.d0 * H) - SIN(H)))
end
;-------------------------------------------------------------------------------
function DP33,H
return,(3.d0 * sqrt(5.d0) / 8.d0 / sqrt(2.d0 ) * (COS(H) - COS(3.d0 * H)))
end
;-------------------------------------------------------------------------------
function DP40,H
return,(-5.d0 / 16.d0 * (7.d0 * SIN(4.d0 * H) + 2.d0 * SIN(2.d0 * H)))
end
;-------------------------------------------------------------------------------
function DP41,H
return,(sqrt(5.d0) / 4.d0 / sqrt(2.d0 ) * (7.d0 * COS(4.d0 * H)+COS(2.d0*H)))
end
;-------------------------------------------------------------------------------
function DP42,H
return,(sqrt(5.d0)/8.d0*(7.d0*SIN(4.d0*H)-2.d0*SIN(2.d0*H)))
end
;-------------------------------------------------------------------------------
function DP43,H
return,(sqrt(35.d0)/4.d0/sqrt(2.d0)*(COS(2.d0*H)-COS(4.d0*H)))
end
;-------------------------------------------------------------------------------
function DP44,H
return,(sqrt(35.d0)/16.d0*(2.d0*SIN(2.d0*H)-SIN(4.d0*H)))
end
;;-------------------------------------------------------------------------------
;function DP50,H
;return,(15.d0/128.d0*(-21.d0*SIN(5.d0*H)-7.d0*SIN(3.d0*H)-2.d0*SIN(H)))
;end
;;-------------------------------------------------------------------------------
;function DP51,H
;return,(sqrt(15.d0)/128.d0*(105.d0*COS(5.d0*H)+21.d0*COS(3.d0*H)+2.d0*COS(H)))
;end
;;-------------------------------------------------------------------------------
;function DP52,H
;return,sqrt(3.d0)*sqrt(35.d0)/64.d0*(15.d0*SIN(5.d0*H)-3.d0*SIN(3.d0*H)-2.d0*SIN(H))
;end
;;-------------------------------------------------------------------------------
;function DP53,H
;return,3.d0*sqrt(35.d0)/128.d0/sqrt(2.d0)*(-15.d0*COS(5.d0*H)+13.d0*COS(3.d0*H)+2.d0 *COS(H))
;end
;;-------------------------------------------------------------------------------
;function DP54,H
;return,(3.d0*sqrt(35.d0)/128.d0*(-5.d0*SIN(5.d0*H)+9.d0*SIN(3.d0*H)-2.d0*SIN(H)))
;end
;;-------------------------------------------------------------------------------
;function DP55,H
;return,(15.d0*sqrt(7.d0)/128.d0/sqrt(2.d0)*(COS(5.d0*H)-3.d0*COS(3.d0* H)+2.d0 *COS(H)))
;end
;;-------------------------------------------------------------------------------
;+
; magnetic_field is a procedure that <behavior desc here>
;
; :Uses:
; coeff, crt_sheet
;
; :Params:
; rr0: in, required, type=float
; rayon
; th0: in, required, type=float
; colatitude
; ph0: in, required, type=float
; azimuth
; b_sph: in, required, type=sometype
; A parameter named b_sph
; b: in, required, type=sometype
; A parameter named b
; br: in, required, type=sometype
; A parameter named br
;
; :Keywords:
; mfl_model: in, optional, type=sometype
; default 'SPV'
; deg: in, optional, type=sometype
; A keyword named deg
; nocrt: in, optional, type=sometype
; A keyword named nocrt
;-
pro MAGNETIC_FIELD, rr0, th0, ph0, b_sph, b, br, mfl_model=mfl_model, deg=deg, nocrt=nocrt
if ~keyword_set(mfl_model) then mfl_model='SPV'
nr = n_elements(rr0)
nt = n_elements(th0)
np = n_elements(ph0)
nk = [nr,nt,np]
nk = nk(sort(nk))
if nk(0)*nk(1)*nk(2) eq 1 then begin
nn = 1
endif else begin
if nk(0)*nk(1) eq 1 then begin
nn = nk(2)
endif else begin
if nk(0) eq 1 then begin
if nk(1) eq nk(2) then begin
nn = nk(1)
endif else message,'dimensions do not agree, aborting.'
endif else begin
if nk(0) eq nk(1) and nk(1) eq nk(2) then begin
nn = nk(0)
endif else message,'dimensions do not agree, aborting.'
endelse
endelse
endelse
if nr eq 1 then r = rebin(reform(rr0,1),nn) else r=rr0
if nt eq 1 then theta = rebin(reform(th0,1),nn) else theta=th0
if np eq 1 then phi = rebin(reform(ph0,1),nn) else phi=ph0
b_sph=dblarr(3,nn)
if keyword_set(deg) then begin
rtheta=double(theta)*!dpi/180.d0
rphi=double(phi)*!dpi/180.d0
endif else begin
rtheta=double(theta)
rphi=double(phi)
endelse
r = double(r)
COEFF,g,h,n,dip,crt,tdip,mfl_model=mfl_model
c=cos( rebin(reform(dindgen(5)+1.,5,1),5,nn) * rebin(reform(rphi,1,nn),5,nn) )
s=sin( rebin(reform(dindgen(5)+1.,5,1),5,nn) * rebin(reform(rphi,1,nn),5,nn) )
b_sph(0,*) = 2.d0 / r^3.d0 * $
( P10(rtheta)* G(1,0) $
+ P11(rtheta)*(G(1,1)*c(0,*)+H(1,1)*s(0,*)) $
) $
+ 3.d0 / r^4.d0 * $
( P20(rtheta)* G(2,0) $
+ P21(rtheta)*(G(2,1)*c(0,*)+H(2,1)*s(0,*)) $
+ P22(rtheta)*(G(2,2)*c(1,*)+H(2,2)*s(1,*)) $
) $
+ 4.d0 / r^5.d0 * $
( P30(rtheta)* G(3,0) $
+ P31(rtheta)*(G(3,1)*c(0,*)+H(3,1)*s(0,*)) $
+ P32(rtheta)*(G(3,2)*c(1,*)+H(3,2)*s(1,*)) $
+ P33(rtheta)*(G(3,3)*c(2,*)+H(3,3)*s(2,*)) $
) $
+ 5.d0 / r^6.d0 * $
( P40(rtheta)* G(4,0) $
+ P41(rtheta)*(G(4,1)*c(0,*)+H(4,1)*s(0,*)) $
+ P42(rtheta)*(G(4,2)*c(1,*)+H(4,2)*s(1,*)) $
+ P43(rtheta)*(G(4,3)*c(2,*)+H(4,3)*s(2,*)) $
+ P44(rtheta)*(G(4,4)*c(3,*)+H(4,4)*s(3,*)) $
); $
; + 6.d0 / r^7.d0 * $
; ( P50(rtheta)* G(5,0) $
; + P51(rtheta)*(G(5,1)*c(0,*)+H(5,1)*s(0,*)) $
; + P52(rtheta)*(G(5,2)*c(1,*)+H(5,2)*s(1,*)) $
; + P53(rtheta)*(G(5,3)*c(2,*)+H(5,3)*s(2,*)) $
; + P54(rtheta)*(G(5,4)*c(3,*)+H(5,4)*s(3,*)) $
; + P55(rtheta)*(G(5,5)*c(4,*)+H(5,5)*s(4,*)) $
; )
b_sph(1,*) = - 1.d0 / r^3.d0 * $
( DP10(rtheta)* G(1,0) $
+ DP11(rtheta)*(G(1,1)*c(0,*)+H(1,1)*s(0,*)) $
) $
- 1.d0 / r^4.d0 * $
( DP20(rtheta)* G(2,0) $
+ DP21(rtheta)*(G(2,1)*c(0,*)+H(2,1)*s(0,*)) $
+ DP22(rtheta)*(G(2,2)*c(1,*)+H(2,2)*s(1,*)) $
) $
- 1.d0 / r^5.d0 * $
( DP30(rtheta)* G(3,0) $
+ DP31(rtheta)*(G(3,1)*c(0,*)+H(3,1)*s(0,*)) $
+ DP32(rtheta)*(G(3,2)*c(1,*)+H(3,2)*s(1,*)) $
+ DP33(rtheta)*(G(3,3)*c(2,*)+H(3,3)*s(2,*)) $
) $
- 1.d0 / r^6.d0 * $
( DP40(rtheta)* G(4,0) $
+ DP41(rtheta)*(G(4,1)*c(0,*)+H(4,1)*s(0,*)) $
+ DP42(rtheta)*(G(4,2)*c(1,*)+H(4,2)*s(1,*)) $
+ DP43(rtheta)*(G(4,3)*c(2,*)+H(4,3)*s(2,*)) $
+ DP44(rtheta)*(G(4,4)*c(3,*)+H(4,4)*s(3,*)) $
); $
; - 1.d0 / r^7.d0 * $
; ( DP50(rtheta)* G(5,0) $
; + DP51(rtheta)*(G(5,1)*c(0,*)+H(5,1)*s(0,*)) $
; + DP52(rtheta)*(G(5,2)*c(1,*)+H(5,2)*s(1,*)) $
; + DP53(rtheta)*(G(5,3)*c(2,*)+H(5,3)*s(2,*)) $
; + DP54(rtheta)*(G(5,4)*c(3,*)+H(5,4)*s(3,*)) $
; + DP55(rtheta)*(G(5,5)*c(4,*)+H(5,5)*s(4,*)) $
; )
b_sph(2,*) = 0
ww = where(rtheta ne 0 and rtheta ne !pi,cnt)
if cnt ne 0 then begin
b_sph(2,ww) = ( 1.d0 / r(ww)^3.d0 * $
P11(rtheta(ww))*(G(1,1)*s(0,ww)-H(1,1)*c(0,ww)) $
+ 1.d0 / r(ww)^4.d0 * $
( P21(rtheta(ww))*(G(2,1)*s(0,ww)-H(2,1)*c(0,ww)) $
+2*P22(rtheta(ww))*(G(2,2)*s(1,ww)-H(2,2)*c(1,ww)) $
) $
+ 1.d0 / r(ww)^5.d0 * $
( P31(rtheta(ww))*(G(3,1)*s(0,ww)-H(3,1)*c(0,ww)) $
+2*P32(rtheta(ww))*(G(3,2)*s(1,ww)-H(3,2)*c(1,ww)) $
+3*P33(rtheta(ww))*(G(3,3)*s(2,ww)-H(3,3)*c(2,ww)) $
) $
+ 1.d0 / r(ww)^6.d0 * $
( P41(rtheta(ww))*(G(4,1)*s(0,ww)-H(4,1)*c(0,ww)) $
+2*P42(rtheta(ww))*(G(4,2)*s(1,ww)-H(4,2)*c(1,ww)) $
+3*P43(rtheta(ww))*(G(4,3)*s(2,ww)-H(4,3)*c(2,ww)) $
+4*P44(rtheta(ww))*(G(4,4)*s(3,ww)-H(4,4)*c(3,ww)) $
) $
; + 1.d0 / r(ww)^7.d0 * $
; ( P51(rtheta(ww))*(G(5,1)*s(0,ww)-H(5,1)*c(0,ww)) $
; +2*P52(rtheta(ww))*(G(5,2)*s(1,ww)-H(5,2)*c(1,ww)) $
; +3*P53(rtheta(ww))*(G(5,3)*s(2,ww)-H(5,3)*c(2,ww)) $
; +4*P54(rtheta(ww))*(G(5,4)*s(3,ww)-H(5,4)*c(3,ww)) $
; +5*P55(rtheta(ww))*(G(5,5)*s(4,ww)-H(5,5)*c(4,ww)) $
; ) $
)/sin(rtheta)
endif
if ~keyword_set(nocrt) then CRT_SHEET,r,rtheta,rphi,crt,tdip,br else br=0.d0
;CRT_SHEET,r,theta,phi,crt,tdip,br
b_sph=b_sph+br*1.d-5
b=sqrt(total(b_sph^2.,1))
return
end