; NAME :	sunparam
; PURPOSE :
;	calculate solar parameters 
; CATEGORY :
; 	solar magnetic field computation package (MAGPACK2)
; CALLING SEQUENCE :
;	sunparam, tsys, date, time, rsun, p, b0, cl0, doy
; INPUTS :
;	tsys = 'GMT+hh', 'JST',,,etc.,  which indicates a timezone
;	date = 'yy/mm/dd', time = 'hh:mm:ss'
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
;	none
; OUTPUTS :
;     rsun=  semidiameter in arcsec
;     p   =  position angle of sun's pole in deg.
;     b0  =  heliographic latitude of disk center in deg.
;     cl0 =  carrington longitude of the meridian in deg.
;     doy =  day of year
; COMMON BLOCKS : 	none
; SIDE EFFECTS : 	none
; RESTRICTIONS :	none
; PROCEDURE :
;
; MODIFICATION HISTORY :
;	originally coded by K.Shibasaki 10/28/83 as koyomi/zs for OAO Melcom 70
;	T.Sakurai, 1992
;	           1993.2 tsys is added to argument list
;
;------------------------------------------------------------------------
;
      pro sunparam, tsys, date, time, rsun, p, b0, cl0, doy

	tntbl=['GMT-12', 'GMT-11', 'HST', 'YST', 'PST', $
               'MST', 'CST', 'EST', 'AST', 'GMT-3',	$
               'GMT-2', 'GMT-1', 'GMT', 'MET', 'EET',	$
               'GMT+3', 'GMT+4', 'GMT+5', 'GMT+6', 'GMT+7',	$
               'WST', 'JST', 'EST', 'GMT+11', 'NZST' ]
;  daylight saving time
	tdtbl=['***-12', '***-11', 'HDT', 'YDT', 'PDT',	$
               'MDT', 'CDT', 'EDT', 'ADT', '***-3',		$
               '***-2', '***-1', 'BST', 'METDST', 'EETDST',	$
               '***+3', '***+4', '***+5', '***+6', '***+7',	$
               '***+8', '***+9', '***+10', '***+11', 'NZDT' ]

	rd=3.14159265/180.

	if n_params() eq 0 then begin
	  tsys=' '
	  print,'enter timezone (GMT+hh, JST, etc.)'
	  read, tsys
	  date = ' '
	  print, 'enter date (yy/mm/dd)'
	  read, date
	  time = ' '
	  print, 'enter time (hh:mm:ss)'
	  read, time
	endif

	tsysu=strtrim(strupcase(tsys),2)

	if tsysu eq 'UT' or tsysu eq 'GMT' then $
		tdif = 0 else begin
	if strmid(tsysu,0,2) eq 'UT' then $
		tdif=fix(strmid(tsysu,2,4)) else begin
	if strmid(tsysu,0,3) eq 'GMT' then $
		tdif=fix(strmid(tsysu,3,3)) else begin
	for tdif=-12, 12 do begin
		if strmid(tsysu,0,strlen(tntbl(tdif+12))) eq tntbl(tdif+12) $
		then goto, found
	endfor
	for tdif=-11, 13 do begin
		if strmid(tsysu,0,strlen(tdtbl(tdif+11))) eq tdtbl(tdif+11) $
		then goto, found
	endfor
	print,'(proc.suncal) unknown timezone:',tsys,', assume UT'
	tdif=0
	endelse
	endelse
	endelse
found:

	year	= fix(strmid(date,0,2))
	month	= fix(strmid(date,3,2))
	day	= fix(strmid(date,6,2))
	hour	= fix(strmid(time,0,2)) - tdif
	min	= fix(strmid(time,3,2))
	sec	= fix(strmid(time,6,2))

	year= (year mod 100)+1900
	ut= ( sec/60.0 + min )/60.0 + hour
	it1=367*(year-1900)
	it2=((month+9)/12+year)*7/4
	it3=275*month/9
	t1=(it1-it2+it3+day) +ut/24.0
	t =(t1+3293.5)/36525.0
	tt=t1+20093.5
	gml=(279.691+36000.769*t)*rd
	g =(358.476+35999.050*t)*rd
	tgl=gml+(1.919-0.0048*t)*rd*sin(g)+0.020*rd*sin(2.0*g)
	rlg=0.00003+(0.000018*t-0.00727)*cos(g)-0.00009*cos(2.0*g)
	dist =10.0^(rlg) ;     dist   distance in au
	rsun =961.18/dist
	e =23.441*rd
	ai =7.25*rd
	omg=(73.739167+0.013958*tt/365.25)*rd
	b0 =asin(sin(tgl-omg)*sin(ai))
	cbclm=-cos(tgl-omg)
	cbslm=-sin(tgl-omg)*cos(ai)
	xlm=atan(cbslm,cbclm)
	cl0=(360.0-360.0/25.38*tt)*rd+xlm
	p =-atan(cos(tgl)*tan(e))-atan(cos(tgl-omg)*tan(ai))
	b0=b0/rd
	cl0=cl0/rd
	lc=fix(cl0/360.0)
	cl0=cl0-360.0*float(lc-1)
	p =p /rd
	idoy1=(month+9)/12*(1+((year mod 4)+2)/3)
	doy =float(it3-idoy1+day-30)+ut/24.0
;
	if n_params() eq 0 then begin
	  print, date,', ',time,' ',tsysu
	  print, 'distance to the sun  =',dist,' a.u.'
	  print, 'solar apparent radius=',rsun,' arcsec'
	  print, 'p,b0,cl0=',p,b0,cl0,' degrees'
	  print, 'doy=',doy
	endif
	return
end