crt_sheet.pro (Kronos IDL library)

;+
; Contains the crt_sheet procedure
;
; :Author:
;	 Baptiste Cecconi
;
; :History:
;	 2008/02/08: verification avec CRT_SHEET de SERPE. -> ok en simple précision
;                    + passage en double précision.
;
;	 2008/02/20: Last Edit
;-
;
;+
; Le but de ce programme est de calculer le champ magnetique du au disque de 
; courant. Pour cela il faut :
;
; 1. Transformer les coordonnees spheriques geographiques en coordonnees 
;    cylindriques magnetiques
;
; 2. Calculer le champ
;
; 3. Transformer le champ donnees en coordonnees cylindriques magnetiques 
;    en champ en coordonnees spheriques geographiques
;
; Les etapes 1 et 3 se font en passant par l'intermediaire de coordonnees 
; cartesiennes
;
; :Params:
;	 r1:	in, required, type=sometype
;		 Coordonnees spheriques geographiques du point dont on desire
;	 t1:	in, required, type=sometype
;		 Coordonnees spheriques geographiques du point dont on desire
;	 p1:	in, required, type=sometype
;		 Coordonnees spheriques geographiques du point dont on desire
;	 crt:  in, required, type=sometype
;		 parametre du courant
;	 t:	 in, required, type=sometype
;		 matrice de changement de base
;	 bres: out, required, type=sometype
;		 Valeurs du champ en ce point
;
; :Keywords:
;	 test: in, optional, type=sometype
;		 A keyword named test
;-
pro CRT_SHEET,r1,t1,p1,crt,t,bres,test=test

nn = n_elements(r1)
tr=dblarr(3,3,nn) 
tr(0,0,*)=sin(double(t1))*cos(double(p1))
tr(1,0,*)=sin(double(t1))*sin(double(p1))
tr(2,0,*)=cos(double(t1))
tr(0,1,*)=cos(double(t1))*cos(double(p1))
tr(1,1,*)=cos(double(t1))*sin(double(p1))
tr(2,1,*)=-sin(double(t1))
tr(0,2,*)=-sin(double(p1))
tr(1,2,*)=cos(double(p1))
tr(2,2,*)=0.

if nn gt 1 then invtr = transpose(tr,[1,0,2]) else invtr = transpose(tr)
; tr is a rotation matrix, so invert==transpose
invt  = transpose(double(t))

ttr = reform(t#reform(tr,3,3l*nn),3,3,nn)  ; ok, checked with old crt_sheet 
rr1 = rebin(reform(double(r1),1,nn),3,nn)
res = reform(ttr(*,0,*)*rr1,3,nn)

r2=sqrt(res(0,*)^2.+res(1,*)^2.)  ; ok, checked with old crt_sheet 

crt = double(crt)

;for i=0l,nn-1l do begin
;  print,r2(i),res(2,i)
;  ;print,ttr(*,*,i)
;endfor

	f1=sqrt((res(2,*)-crt(2))^2+crt(1)^2)
	f2=sqrt((res(2,*)+crt(2))^2+crt(1)^2)
	brp2=r2/2.*(1./f1-1./f2)
	bzp2=2*crt(2)/sqrt(res(2,*)^2.+crt(1)^2.)- $
	r2^2/4.*((res(2,*)-crt(2))/f1^3.-(res(2,*)+crt(2))/f2^3.)

br2 = brp2*0.d0
bz2 = bzp2*0.d0

if keyword_set(test) then print,crt(0)
if keyword_set(test) then print,r2
if keyword_set(test) then print,crt(2)
if keyword_set(test) then print,abs(res(2,*))

w1 = where(r2 lt crt(0),cnt1)
if cnt1 ne 0 then begin
if keyword_set(test) then print,"case 1",w1
	f1(w1)=sqrt((res(2,w1)-crt(2))^2+crt(0)^2)
	f2(w1)=sqrt((res(2,w1)+crt(2))^2+crt(0)^2)
	br2(w1)=crt(3)*(r2(w1)/2.*(1./f1(w1)-1./f2(w1))-brp2(w1))
	bz2(w1)=crt(3)*(2*crt(2)/sqrt(res(2,w1)^2.+crt(0)^2.)- $
			 r2(w1)^2/4.*((res(2,w1)-crt(2))/f1(w1)^3.-(res(2,w1)+ $
			 crt(2))/f2(w1)^3.)-bzp2(w1))
endif

w2 = where(r2 ge crt(0) and abs(res(2,*)) gt crt(2), cnt2)
if cnt2 ne 0 then begin
if keyword_set(test) then print,"case 2",w2
	f1(w2)=sqrt((res(2,w2)-crt(2))^2+r2(w2)^2)
	f2(w2)=sqrt((res(2,w2)+crt(2))^2+r2(w2)^2)
	br2(w2)=crt(3)*(((f1(w2)-f2(w2)+abs(res(2,w2))/res(2,w2)*2*crt(2))/r2(w2)-$ 
			 crt(0)^2*r2(w2)/4.*(1./f1(w2)^3-1./f2(w2)^3))-brp2(w2))
	bz2(w2)=crt(3)*(2*crt(2)/sqrt(res(2,w2)^2.+r2(w2)^2.)- $
			 crt(0)^2/4.*((res(2,w2)-crt(2))/f1(w2)^3.-(res(2,w2)+crt(2))/f2(w2)^3.)-bzp2(w2))
endif

w3 = where(r2 ge crt(0) and abs(res(2,*)) le crt(2), cnt3)
if cnt3 ne 0 then begin
if keyword_set(test) then print,"case 3",w3
	f1(w3)=sqrt((res(2,w3)-crt(2))^2+r2(w3)^2)
	f2(w3)=sqrt((res(2,w3)+crt(2))^2+r2(w3)^2)
	br2(w3)=crt(3)*(((f1(w3)-f2(w3)+2.*res(2,w3))/r2(w3)-crt(0)^2*r2(w3)/4.*(1./f1(w3)^3-1./f2(w3)^3))-brp2(w3))
	bz2(w3)=crt(3)*(2*crt(2)/sqrt(res(2,w3)^2.+r2(w3)^2.)- $
	crt(0)^2/4.*((res(2,w3)-crt(2))/f1(w3)^3.-(res(2,w3)+crt(2))/f2(w3)^3.)-bzp2(w3))
endif

; Br2 and Bz2  : ok, checked with old crt_sheet 

bres = [br2*res(0,*)/r2,br2*res(1,*)/r2,bz2]


for i=0l,nn-1l do bres(*,i) = invtr(*,*,i)#invt#bres(*,i)
return
end