;+ ; 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' or 'yyyymmdd', 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 ; 1999.1 format of year can be 'yyyymmdd' ;- ;------------------------------------------------------------------------ ; 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.14159265358979d0/180.0d0 if n_params() eq 0 then begin tsys=' ' print,'enter timezone (GMT+hh, JST, etc.)' read, tsys date = ' ' print, 'enter date (yy/mm/dd or yyyymmdd)' 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: if (strmid(date,2,1) eq '/') then begin year = fix(strmid(date,0,2)) if (year le 50) then year=year+2000 else year=year+1900 month = fix(strmid(date,3,2)) day = fix(strmid(date,6,2)) endif else begin year = fix(strmid(date,0,4)) month = fix(strmid(date,4,2)) day = fix(strmid(date,6,2)) endelse hour = double(strmid(time,0,2)) - tdif min = double(strmid(time,3,2)) sec = double(strmid(time,6,2)) ut= ( sec/60.0d0 + min )/60.0d0 + hour it1=367L*(year-1900) it2=((month+9)/12+year)*7/4 it3=275L*month/9 t1=double(it1-it2+it3+day) +ut/24.0d0 t =(t1+3293.5d0)/36525.0d0 tt=t1+20093.5d0 gml=(279.691d0 + 36000.769d0*t)*rd g =(358.476d0 + 35999.050d0*t)*rd tgl=gml+(1.919d0-0.0048d0*t)*rd*sin(g)+0.020d0*rd*sin(2.0d0*g) rlg=0.00003d0+(0.000018d0*t-0.00727d0)*cos(g)-0.00009d0*cos(2.0d0*g) dist = float(10.0d0^(rlg)) ; dist distance in au rsun = float(961.18d0/dist) e =23.441d0*rd ai =7.250d0*rd omg=(73.739167d0+0.013958d0*tt/365.25d0)*rd b0 =asin(sin(tgl-omg)*sin(ai)) cbclm=-cos(tgl-omg) cbslm=-sin(tgl-omg)*cos(ai) xlm=atan(cbslm,cbclm) cl0=(360.0d0-360.0d0/25.38d0*tt)*rd+xlm p =-atan(cos(tgl)*tan(e))-atan(cos(tgl-omg)*tan(ai)) b0=float(b0/rd) cl0=cl0/rd lc=fix(cl0/360.0d0) cl0=float(cl0-360.0d0*(lc-1)) p = float(p/rd) idoy1=(month+9)/12*(1+((year mod 4)+2)/3) doy =float(it3-idoy1+day-30 +ut/24.0d0) ; 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