;+
; PURPOSE : calculate Fractional Carrington longitude (and B0, Rsun, P angle;           corresponding to CMP time (UT)
; CALLING SEQUENCE :
;       cl0_yyyymmddtime,year,month,day,hrs,cl0
; INPUTS : year,month,day,hrs (in fractinal hours)
; OUTPUTS : cl0 =  carrington longitude of the CMP
; MODIFICATION HISTORY :
;       First version: 20090531 by X. Zhao. The formulae used are the 
;                      same as what in sunparams.pro by T.Sakurai, 2000.1
;-
;------------------------------------------------------------------------
;
pro cl0_date,date,cl0,PLT=plt

    if n_params( ) lt 1 then begin
       print,'cl0_date,date,cl0,/plt'
       return
    endif
    year_sec_date,date,year,month,day,hh,mm,ss
    hrs=hh+mm/60.+ss/3600.
    rd=3.14159265358979d0/180.0d0
    it1=367L*(year-1900)
    it2=((month+9)/12+year)*7/4
    it3=275L*month/9
    t1=double(it1-it2+it3+day) +hrs/24.0d0
; t1 = JD -2418313.5
    t =(t1+3293.5d0)/36525.0d0
; t = (JD -2415020)/36525
    tt=t1+20093.5d0
; tt = JD -2398220
    gml=(279.691d0 + 36000.769d0*t)*rd
    g =(358.476d0 + 35999.050d0*t)*rd
    tgl=gml+(1.919d0-0.0048d0*t)*rd*sin(g)+0.020d0*rd*sin(2.0d0*g)
    rlg=0.00003d0+(0.000018d0*t-0.00727d0)*cos(g)-0.00009d0*cos(2.0d0*g)
    dist = float(10.0d0^(rlg)) ;     dist   distance in au
    rsun = float(961.18d0/dist)
    e =23.441d0*rd
    ai =7.250d0*rd
    omg=(73.739167d0+0.013958d0*tt/365.25d0)*rd
    b0 =asin(sin(tgl-omg)*sin(ai))
    b0=float(b0/rd)
    cbclm=-cos(tgl-omg)
    cbslm=-sin(tgl-omg)*cos(ai)
    xlm=atan(cbslm,cbclm)
    cl0=(360.0d0-360.0d0/25.38d0*tt)*rd+xlm
    cl0=cl0/rd
    lc=fix(cl0/360.0d0)
    cl0=float(cl0-360.0d0*(lc-1))
    p =-atan(cos(tgl)*tan(e))-atan(cos(tgl-omg)*tan(ai))
    p = float(p/rd)
    idoy1=(month+9)/12*(1+((year mod 4)+2)/3)
    doy =float(it3-idoy1+day-30 +hrs/24.0d0)
    if keyword_set(plt) then begin
       date_year_sec,year,month,day,hh,mm,ss,date
       print,'date:',date
       print,'cl0,rsun,b0,p:',cl0,rsun,b0,p
    endif

end