pro suneph,Date,Time,SUN,path=path, universal=universal

;+
; This procedure reads contents of the file containing solar data
; (e.g. for year 1993 - file "sol93.dat") and interpolates data to the
; observation time. File sol**.dat is assumed to be placed into the
; directory described by the environment variable "astr_data". If it is
; not the case, you must specify keyword parameter "path".
;
; Input arguments Date and Time are strings,
;    e.g. Date=12 07 94 - July 12, 1994;
;         Time=21 23 45.385 - 21 hour 23 min 45.385 sec UT.
;    Second fractions are not obliged.
;
; All the output values are collected into the structure SUN
; named SOL_EPHEMERIDE.
;
; Output values are concerned to the observation time:
;       SUN.H	 - hour angle;
;       SUN.Decl - declination;
;	SUN.W0	 - angular velocity of the diurnal rotation of
;                the Earth, rad/sec;
;	SUN.Tcul - culmination time for the given day in hours;
;	SUN.R	 - optical radius of the Sun;
;	SUN.Dp	 - position angle, i.e the angle between the
; 		 polar axes of the Sun and the meridium being
;		 counted to the east from the meridium (diurnal parallel);
;	SUN.B0   - altitude of the solar center;
;	SUN.Karr - Karrington longitude of the solar center.
; All the output values are floating-point, double precision and are
; measured (except for Tcul and W0) in radians.
;
;	SUN.Current_date - buffer strings containing summary information read
; from file (about three consequtive days around current date).
; SUN.Current_date(3) contains current year.
;
;-

DDTOR=!DPi/180

sz=size(SUN)
sun_type=sz(n_elements(sz)-2)

if sun_type eq 8 then if strupcase(tag_names(SUN,/str)) eq $
	'SOL_EPHEMERIDE' then begin

SUN.Date=Date
SUN.Time=Time
Month=strmid(SUN.Current_date(1),5,2)
Nday=strmid(SUN.Current_date(1),7,2)
if strmid(Nday,0,1) eq ' ' then Nd='0'+strmid(Nday,1,1) else Nd=Nday
if strmid(Month,0,1) eq ' ' then Mn='0'+strmid(Month,1,1) else Mn=Month
;Date_Buf=Nd+' '+Mn+' '+strmid(SUN.Current_date(3),2,2)
Date_Buf=Nd+'/'+Mn+'/'+strmid(SUN.Current_date(3),2,2)
Date_ = Date

strput, Date_, '/', 2
strput, Date_, '/', 5


	If Date_ eq Date_Buf then begin
Current_date=SUN.Current_date
goto, Current
	endif

		endif

New:
		Year=strmid(Date,6,2)

Current_date=['','','',Year]

if Year gt 50 then Current_date(3)='19'+Year else Current_date(3)='20'+Year
if n_elements(path) le 0 then path=getenv('astr_data')

	CASE !version.OS OF

'windows':	Delim='\'

'Win32':	Delim='\'

ELSE:		Delim='/'

	ENDCASE



Name=path+Delim+'sol'+Year+'.dat'
NdayI=fix(strmid(Date,0,2))
MonthI=fix(strmid(Date,3,2))

if(float(Year)/4 eq fix(Year)/4 and Year ne '00') then Mvys=1 else Mvys=0
		M=MonthI
IF(NdayI EQ 1) THEN BEGIN MonthP=MonthI-1
    IF(M EQ 5 OR M EQ 7 OR M EQ 10 OR M EQ 12)THEN NdayP=30 $
    ELSE IF (M EQ 3 AND Mvys EQ 1) THEN NdayP=29 $
    ELSE IF (M EQ 3 AND Mvys EQ 0) THEN NdayP=28 ELSE NdayP=31
ENDIF ELSE BEGIN NdayP=NdayI-1
		 MonthP=MonthI
ENDELSE
print,'reading the file ',Name

if !d.flags and 2l^16 ne 0 then WIDGET_CONTROL,/hour

	openr,lun,Name,ERROR=err,/get_lun
	IF(err NE 0) then PRINTF,-2,!ERR_STRING

header=strarr(3)
status=fstat(lun)
	readf,lun,header
	point_lun,-lun,header_length

date0='01 01 '+year
t_dif=time_difference(date0,time,date,time)/3600d0/24

	IF t_dif GT 10 THEN BEGIN

Last_record=status.size-200

	point_lun,lun,	Last_record

a=bytarr(200)

	readu,lun,a

ind=where(a eq '0a'xb)

	point_lun,lun, Last_record+ind(0)+1

ws=''

	readf,lun,ws
	point_lun,-lun, Last_pointer

Last_Date=strmid(ws,7,2)+' '+strmid(ws,5,2)+' '+Year

	if time_difference(date,time,Last_Date,time) lt 0 then begin
	free_lun,lun
	message,'Error! Ephemeride for this date is not available.'
	endif

ndays=fix(time_difference(date0,time,Last_Date,time)/3600d0/24+0.5)+1

string_length=(Last_pointer-header_length)/float(ndays+1)

first_step= long((t_dif-5)*string_length)+header_length

	point_lun,lun,	first_step

temp=bytarr(150)

	readu,lun,temp

ind=where(temp eq '0a'xb)

	point_lun,lun, first_step+ind(0)+1

	ENDIF

Data_string=''

			I=0

		WHILE I lt 3 DO BEGIN
	readf,lun,Data_string
	Current_date(I)=Data_string
		Month=fix(strmid(Current_date(0),5,2))
		Nday=fix(strmid(Current_date(0),7,2))
IF(Month EQ MonthP AND Nday EQ NdayP OR MonthP EQ 0) THEN I=I+1
		ENDWHILE

	free_lun,lun

Current:
	Declg=strmid(Current_date,12,2)
	DSign=strmid(Current_date,11,1)
	Declm=strmid(Current_date,15,2)
	Decls=strmid(Current_date,18,4)
	Hculm=strmid(Current_date,24,2)
	Mculm=strmid(Current_date,27,2)
	Sculm=strmid(Current_date,30,4)
	Rsun=strmid(Current_date,35,5)*DDTOR/60
	Dpar=strmid(Current_date,42,5)*DDTOR
	B0in=strmid(Current_date,49,4)*DDTOR
	Karr=strmid(Current_date,55,6)*DDTOR
	Ddelta=strmid(Current_date,62,6)*DDTOR/3600d0*24d0


Tculm=HMS(Hculm,Mculm,Sculm)
	Tcul=Tculm(1)
Delta=HMS(Declg,Declm,Decls)*DDTOR
for I=0,2 do IF(Dsign(I) EQ '-') then Delta(I)=-Delta(I)

	Tobs=HMS(Time)

	Tobs=Tobs-24*(Tobs ge (Tcul+12))*(1-keyword_set(universal))

IF((Tobs-Tcul) LE 0) THEN Tdiur=Tculm(1)-Tculm(0)+24 ELSE $
	Tdiur=Tculm(2)-Tculm(1)+24
		W0=2*!DPi/Tdiur
		Hangle=W0*(Tobs-Tcul)
		W0=W0/3600d0

Direction=Hangle/(2*!DPi)

	IF(Direction LE 0) THEN BEGIN
Decl=Delta(1)+Direction*Ddelta(1)-Direction^2*(DDelta(1)-DDelta(0))/2
Dp=Dpar(1)+Direction*(Dpar(1)-Dpar(0))
B0=B0in(1)+Direction*(B0in(1)-B0in(0))
;Klong=Karr(1)+Direction*(Karr(1)-Karr(0))
Klong=Karr(1)+Direction*(Karr(1)-Karr(0)-2*!DPi*(Karr(1) gt Karr(0)))
Rsol=(Rsun(1)+Direction*(Rsun(1)-Rsun(0)))
	ENDIF ELSE BEGIN
Decl=Delta(1)+Direction*Ddelta(1)+Direction^2*(DDelta(2)-DDelta(1))/2
Dp=Dpar(1)+Direction*(Dpar(2)-Dpar(1))
B0=B0in(1)+Direction*(B0in(2)-B0in(1))
;Klong=Karr(1)+Direction*(Karr(2)-Karr(1))
Klong=Karr(1)+Direction*(Karr(2)-Karr(1)-2*!DPi*(Karr(2) gt Karr(1)))
Rsol=(Rsun(1)+Direction*(Rsun(2)-Rsun(1)))
	ENDELSE
Klong=Klong+2*!DPi*(Klong lt 0)-2*!DPi*(Klong gt 2*!DPi)

SUN={SOL_EPHEMERIDE, Current_date:Current_date, Date:Date, Time:Time, H:Hangle, $
Decl:Decl, W0:W0, Tcul:Tcul, R:Rsol, Dp:Dp, B0:B0, Karr:Klong}

exit:

end