;
; PURPOSE: calculate theta, phi angles and their functions
; INPUT:   iph,jth
; OUTPUT:  phd,thd,lad,cth,sth,csc,scs  
; WRITTEN: 16Feb1999 Xuepu Zhao
; MODIFIED: 05Dec1999   add RUP (starting from Upper-Right)
; MODIFIED: 01APR2004 	Xuepu Zhao for adding /dbl
;
PRO pzgrid,iph,jth,phd,thd,lad,cth,sth,RUP=rup,DBL=dbl

    if N_params() LT 1 then begin
       print,'PRO pzgrid,iph,jth,phd,thd,lad,cth,sth,/rup,/dbl'
       return
    endif
; calculate cth,sth,phd,thd,lad
    if keyword_set(dbl) then  begin
       phd=dblarr(iph)
       thd=dblarr(jth)
       dph=double(360./iph)
       dcth=double(2.0/jth)
    endif else begin
       phd = fltarr(iph)
       thd = fltarr(jth)
       dph=360./iph
       dcth=2.0/jth
    endelse
    lad = thd
    sth = thd
    cth = thd
; Starting from lower left to upper right
    for i=0,iph-1 do phd(i) = (i+1.0)*dph  ; from left to right
    for j=0,jth-1 do begin
        cth(j)=(j+0.5)*dcth-1.0          ; from south to north but not -1=>1
        thr=acos(cth(j))                 ; because poles are geometrical
        sth(j) = SIN(thr)                ; singular points
        thd(j) = thr/!dtor
        lad(j) = 90-thd(j)
    endfor
; Starting from upper right to lower left
     if keyword_set(rup) then begin
        lphd=phd & lthd=thd & llad=lad & lcth=cth & lsth=sth 
        for i=0,iph-1 do phd(i)=lphd(iph-1-i)
        for j=0, jth-1 do begin
            thd(j)=lthd(jth-1-j)
            lad(j)=llad(jth-1-j)
            cth(j)=lcth(jth-1-j)
            sth(j)=lsth(jth-1-j)
        endfor
     endif

END