; NAME :	stransit
; PURPOSE :
;	calculate the meridian transit time in UT of the sun at a specified
;	date and at a specified longitude
; CATEGORY :
; 	practical astronomy
; CALLING SEQUENCE :
;	tmt = stransit(date, longitude=longitude)
; INPUTS :
;	date given as 'yy/mm/dd'
; OPTIONAL INPUT PARAMETERS : 
;	longitude :	if not specified, assume longitude=0.0
; KEYWORD PARAMETERS :
; OUTPUTS :
;	tmt = meridian transit time of the sun in UT hour
; COMMON BLOCKS : 	none
; SIDE EFFECTS :
; RESTRICTIONS :
; PROCEDURE :
;
; MODIFICATION HISTORY :
;	based on Almanac for Computors (U.S.Naval Observatory)
;	originally coded by K.Shibasaki 10/28/83 as koyomi/zs for OAO Melcom 70
;	T.Sakurai, 1993
;
function stransit, date, longitude=longitude

	rd=3.14159265/180.

	if n_params() eq 0 then begin
	  date = ' '
	  print, 'enter date (yy/mm/dd)'
	  read, date
	  longitude=0.0
	  print, 'enter longitude'
	  read, longitude
	endif

	year	= fix(strmid(date,0,2))
	month	= fix(strmid(date,3,2))
	day	= fix(strmid(date,6,2))

	year= (year mod 100)+1900
	it3=275*month/9
	idoy1=(month+9)/12*(1+((year mod 4)+2)/3)
	doy =float(it3-idoy1+day-30)

	case (year) of
      1982: begin
		x1=3.307
		x2=9.332
		x3=37.329
		x4=0.91747
		end
      1983: begin
		x1=3.562
		x2=9.094
		x3=36.377
		x4=0.91747
		end
      1984: begin
		x1=3.819
		x2=8.855
		x3=35.421
		x4=0.91747
		end
      1985: begin
		x1=3.508
		x2=9.598
		x3=38.394
		x4=0.91747
		end
      1986: begin
		x1=3.761
		x2=9.363
		x3=37.454
		x4=0.91747
		end
      1987: begin
		x1=4.020
		x2=9.122
		x3=36.486
		x4=0.91747
		end
      1988: begin
		x1=4.274
		x2=8.887
		x3=35.546
		x4=0.91747
		end
      1989: begin
		x1=3.539
		x2=9.639
		x3=38.557
		x4=0.91747
		end
	1990: begin
		x1=3.798
		x2=9.397
		x3=37.589
		x4=0.91747
		end
      1991: begin
		x1=4.052
		x2=9.161
		x3=36.645
		x4=0.91748
		end
      1992: begin
		x2=8.92889
		x3=35.7187
		x4=0.917505
		end
      1993: begin
		x1=3.49347
		x2=9.67955
		x3=38.71657
		x4=0.917486
		end
      1994: begin
		doy=doy+365.
		x1=3.49347
		x2=9.67955
		x3=38.71657
		x4=0.917486
		end
	else: begin
		print, '(func.stransit) no ephemeris parameters available',year
		return, 0
		end
	endcase

      th1=1.915*sin((0.9856*doy-x1)*rd)
      th2=0.020*sin((1.9712*doy-2.0*x1)*rd)
      if(year le 1984) then $
        th3=0. $
      else $
        th3=.014*cos((.9856*doy-x1)*rd)

      th =(x2+0.98561*doy+th1+th2+th3)*rd
      atth=atan(tan(th)/x4)/rd
	while(1) do begin
		sa1=abs(atth-th/rd)
		sa2=abs(atth+180.-th/rd)
		if(sa1 lt sa2) then goto, escape
		atth=atth+180.
	endwhile
escape:
;  eqt=equation of time
      eqt=(x3+3.94244*doy-4.*atth)/60.
      tmt=12.-eqt-longitude/15.0
	tmt1=tmt
	if(tmt lt 0) then tmt1=tmt1+24.0
	itmt=fix(tmt1)
	min=fix(60*(tmt1-itmt))
	sec=fix(3600*(tmt1-itmt)-60*min)

	if n_params() eq 0 then begin
	  print, '(func.stransit) date=',date,'  longitude=',longitude
	  print, 'meridian transit time of the sun = ', $
		tmt,' = ',string(format='(i2.2,":",i2.2,":",i2.2)', $
		itmt,min,sec), ' UT'
	endif
	return, tmt
      end