;********************************************************************
;NAME       : polangle (procedure)
;FUNCTION   : calculate the helioparameters P, L0, B0
;DATE       : 91/09/27       1991,7,19
;PROGRAMMER : k.i. 	from fang's prog
;********************************************************************
;The calculation is according to the "Astronomical Formulae for calculators"
;Jean Meeus , Published by Willmann-Bell,Inc.
pro polangle,date,time,P,B0,L0
; date : 'yy/mm/dd'
; time : 'hh:mm:ss' in JST
; P    : polar angle (output)
; B0   : 
; L0   :

;---calculate the Julian Date  , valid for after year 1582
hh=strmid(time,0,2)
mi=strmid(time,3,2)
ss=strmid(time,6,2)
t=(double(hh)+double(mi)/60+double(ss)/3600-9)/24 ; UT

yy=double('19'+strmid(date,0,2))
mo=double(strmid(date,3,2))
dd=double(strmid(date,6,2)) +t


if (mo eq 1) or (mo eq 2) then begin 
	y=yy-1
	m=mo+12
endif else begin
	y=yy
	m=mo
endelse 


a=fix(y/100)
b=2-a+fix(a/4)

JD =long(365.25*y)+long(30.6001*(m+1))+ dd + 1720994.5 + b  ;in Gregotin calendar.

T = (JD - 2415020.0)/36525
deltaT = 0.41 + 1.2053 + 0.4992*T*T  ; deltaT=ET-UT   in minutes
JD=JD+deltaT/60/24

Theta    = (JD-2398220)*360/25.38 mod 360   ; sidereal time 
I     	 = 7.25
K   	 = 74.3646 + 1.395833*T

L        = (279.69668 + 36000.76892*T + 0.0003025*T*T ) mod 360
M        = (358.47583 + 35999.04975*T - 0.000150*T*T ) mod 360 ; sun's mean anomaly
C        = ( (1.919460 - 0.004789*T -0.000014*T*T)*sin(M) $
	   +(0.020094 - 0.000100*T)*sin(2*M)             $
	   +0.000293*sin(3*M) ) mod 360
GLS      = (L+C ) mod 360 ; geocentric longitude of the sun
Omega    = (259.18 - 1934.142*T ) mod 360
Epsilon  = 23.452294-0.0130125*T

Pi       = 3.1415927
Lambda   = (GLS -0.00569 ) mod 360
Lambda1  = (Lambda - 0.00479*sin(Omega*Pi/180) ) mod 360 
x        = atan(-cos(Lambda1*Pi/180)*tan(Epsilon*Pi/180))
y        = atan(-cos((Lambda-K)*Pi/180)*tan(I*Pi/180))
P        = x+y
P        = P*180/Pi   ; in degree
B0       = (asin(sin((Lambda-K)*Pi/180)*sin(I*Pi/180)))*180/Pi mod 360

Eta      = (atan(tan((Lambda-K)*Pi/180)*cos(I*Pi/180)))*180/Pi
sign     = abs(Eta-(Lambda-K-180))
sign1    = abs(Eta-(Lambda-K+180))

if ( (sign gt 90) or (sign1 gt 90)) then $
	Eta=Eta+180
L0       = (Eta-Theta) mod 360
if L0 lt 0 then L0=L0+360     ; L0  0---360  

;print,'P = ',P
;print,'B0= ',B0
;print,'L0= ',L0
return
end