;
; PURPOSE: calculate theta, phi angles and their functions
; INPUT:   jph and ith
; OUTPUT:  cth,sth,thd,phd   0->jph-1, 0->ith-1
; WRITTEN: 16Feb1999 Xuepu Zhao
; MODIFIED: 05Dec1999   add RUP (starting from Upper-Right)
; MODIFIED: 19Apr2003 	Xuepu Zhao for the first value of phd
; MODIFIED: 01APR2004 	Xuepu Zhao for adding /dbl
;
PRO zgrid,jph,ith,phd,thd,lad,cth,sth,csc,scs,CMP=cmp,RUP=rup,DBL=dbl

    if N_params() LT 1 then begin
       print,'PRO zgrid,jph,ith,phd,thd,lad,cth,sth,cmp,scs,/rup,/dbl'
       print,'/rup: thd(i)=lthd(ith-i) & phd(j)=lphd(jph-j)' 
       return
    endif
; calculate cth,sth,thd,phd
    if keyword_set(dbl) then  begin
       thd=dblarr(ith)
       phd=dblarr(jph)
       dcth=double(2.0/ith)
       dph=double(360./jph)
       csc=dblarr(ith)
       cpc=dblarr(jph)
    endif else begin
       thd = fltarr(ith)
       phd = fltarr(jph)
       dcth=2.0/ith
       dph=360./jph
       csc=fltarr(ith)
       cpc=fltarr(jph)
    endelse
    lad = thd
    sth = thd
    cth = thd
    for i=0,ith-1 do begin
        cth(i)=(i+0.5)*dcth-1.0          ; from south to north
        thr=acos(cth(i))
        sth(i) = SIN(thr)
        thd(i) = thr/!dtor
        lad(i) = 90-thd(i)
    endfor
    for j=0,jph-1 do phd(j) = j*dph+dph  ; from left to right
; calculate csc,scs
    if keyword_set(cmp) eq 1 then begin
       scs=cpc
       for i=0,ith-1 do csc(i)=1/sth(i)
       for j=0,jph-1 do begin
           dphr=(phd(j)-cmp)*!dtor
           cpc(j)=cos(dphr)
           scs(j)=1/cpc(j)
       endfor
     endif
     if keyword_set(rup) then begin
        lphd=phd & lthd=thd & llad=lad & lcth=cth & lsth=sth & lscs=scs
        for i=0, ith-1 do begin
            thd(i)=lthd(ith-i)
            lad(i)=llad(ith-1)
            cth(i)=lcth(ith-i)
            sth(i)=lsth(ith-1)
        endfor
        for j=0,jph-1 do begin
            phd(j)=lphd(jph-j)
            scs(j)=lscs(jph-j)
        endfor
     endif

END