;
; WRITTEN: 08JAN2009	Xuepu
;
pro lat4carringtonrate,ix,jy,aa,bb,cc,lat2,ind,SNOG=snog,PEVT=pevt,MEUN=meun,PHIL=phil

    if n_params( ) lt 1 then begin
       print,'lat4carringtonrate,ix,jy,aa,bb,cc,lat2,ind,/snog,/meun,/pevt'
       print,'lat4carringtonrate,ix,jy,aa,bb,cc,lat2,ind'
       print,'Say ix=360 & jy=180. By setting keyparam, aa, bb, cc are given'
       return
    endif
    if keyword_set(meun) then begin   ; Meunier 
       aa=13.562 & bb=-2.04 & cc=-1.4875
    endif
    if keyword_set(phil) then begin      ; Phil, ApJ 241:811-819, 1980
    ;  lim=0.006 
       aa=2.917 & bb=-0.40 & cc=-0.40       ;in mu rad / s, sidereal 
       dd=0.202006                       ;1 deg / day=0.202006 mu rad / s 
       aa=aa/dd & bb=bb/dd & cc=cc/dd        ;in deg / day
       ee=0.930505                       ;synodic/sidereal=0.930505
       aa=aa*ee & bb=bb*ee & cc=cc*ee    ;in synodic
    endif
    if keyword_set(snog) then begin
    ;  lim=0.006 
       aa=2.897 & bb=-0.339 & cc=-0.485      ;in mu rad / s, sidereal 
       dd=0.202006                       ;1 deg / day=0.202006 mu rad / s 
       aa=aa/dd & bb=bb/dd & cc=cc/dd        ;in deg / day
       ee=0.930505                       ;synodic/sidereal=0.930505
       aa=aa*ee & bb=bb*ee & cc=cc*ee    ;in synodic
    endif
    if keyword_set(pevt) then begin
    ;  lim=0.009
       aa=13.51 & bb=-1.72 & cc=-2.31
    endif
  ; print,'aa,bb,cc:',aa,bb,cc
    csr=13.1988                           ;CR rate in deg / day, synodic
    zgrid,ix,jy,phd,thd,lad,cth,sth
  ; slar2=fltarr(jy) & slar4=slar2
    clar2=fltarr(jy) & clar4=clar2
    rrt=clar2
    for jj=0,jy-1 do begin
    ;   slarj=SIN(lad(jj)*!DTOR)
    ;   slar2(jj)=slarj*slarj
    ;   slar4(jj)=slarj^4.
    ;   rrt(jj)=aa+bb*slar2(jj)+cc*slar4(jj)
        clarj=cth(jj)
        clar2(jj)=clarj*clarj
        clar4(jj)=clarj^4.
        rrt(jj)=aa+bb*clar2(jj)+cc*clar4(jj)
    endfor
;   print,rrt                            ; find out lim
;   plot,lad,rrt,xrange=[-90,90],xstyle=1
;   oplot,[-90,90],[csr,csr]
    drrt=ABS(rrt-csr)
    min=MIN(drrt)
    ind=where(drrt eq min,c0)
    if c0 gt 0 then begin
       lat2=thd(ind)-90
    ;  print,'ind:',ind
    ;  print,'lad(ind):',lat2
    ;  print,'rrt(ind): ',rrt(ind)
    endif else print,'No match'

end