PRO zgrid,jph,ith,phd,thd,lad,cth,sth,cmp,scs,RUP=rup

;+
; NAME:
;	zgrid
; PURPOSE:
;	calculate theta, phi angles and their functions
; INPUT:
;	jph 	# longitudes
;			the phi grid is an equally spaced open grid covering
;			(0,360) degrees
;	ith 	# latitudes
;			the theta grid is an equally spaced open grid in
;			cos(theta) covering (-1,1)
;	cmp 	?? some angle in degrees
; OPTIONAL INPUTS:
;	/rup	??
; OUTPUT:
;	phd 	phi   (degrees)
;	thd 	theta (degrees)
;	lad 	latitude (degrees) (=90-thd)
;	cth 	cos(theta)
;	sth 	sin(theta)
; REMARKS:
;	PPH: /rup shouldn't work: the do-loops access arrays
;		with index out of range.
; MODIFICATION HISTORY:
;	16Feb1999 Xuepu Zhao
;	05Dec1999 add RUP (starting from Lower_Left)
;-

IF N_params() LT 2 THEN $
	message, 'SYNTAX - PRO zgrid,jph,ith,phd,thd,lad,cth,sth,cmp,scs,/rup'

; calculate cth,sth,thd,phd

dcth=2.0/ith
thd = fltarr(ith)
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

;cth = gridgen(ith,range=[-1,1],/open)
;thr = acos(cth)
;thd = thr/!dtor
;sth = sin(thr)
;lad = 90-thd

phd = fltarr(jph)
dph = 360./jph
for j=0,jph-1 do phd(j) = j*dph+dph/2.  ; from left to right

;phd = gridgen(jph,range=[0,360],/open)

; calculate csc,scs

if n_elements(cmp) eq 1 then begin
	csc=1/sth
	cpc=cos((phd-cmp)*!dtor)
	scs=1/cpc
	lscs=scs
	if keyword_set(rup) then scs=reverse(scs)
endif

if keyword_set(rup) then begin
	lphd=phd
	lthd=thd
	llad=lad
	lcth=cth
	lsth=sth
	phd=reverse(lphd)
	thd=reverse(lthd)
	lad=reverse(llad)
	cth=reverse(lcth)
	sth=reverse(lsth)
endif

RETURN  &  END