program culplot
use f90_unix_env    ! F compiler only
use f90_kind
use a137
implicit none
!============================================================================
! Site specific:
character(len=*), parameter :: TACTICAL = "c:\tactical\"
character(len=*), parameter :: WGETPATH = "c:\wget\wget"
! OS specific: for Linux, change black to white!
character(len=*),dimension(8),parameter :: shade = (/"black  ", "red    ", &
	"green  ", "blue   ", "cyan   ", "magenta", "yellow ", "orange "/)
! No need to change anything below this line ================================

!integer, parameter :: double = 8, single = 4 ! decomment for g95
! GLOBALS:
integer, parameter :: NP = 100, delay = 90      ! second delay between frames
logical, parameter :: WEB = .false.  ! .true. for cgi, .false. for standalone
real, parameter :: hr=3600.0, r2d=57.29578, lat=-30.23797, glat=-30.23797, &
	altobs=2710.0, long=-70.73369/15.0
real(kind=double) :: FJ
integer :: nC,epoch,yr,mon,day,object,nO,tick
real :: seeingAvg,utnow = 0,ratop,dectop,sunset,sunrise,tucrepv,tucrepm,ts0
real, dimension(1000) :: utC,raC,decC,seeing,elevC
integer, dimension(1000) :: inst
character(len=22) :: cap
character(len=100) :: pastfile, datafile
character(len=16),dimension(20) :: nomO
character(len=44),dimension(20) :: msgO
real,dimension(20) :: raO,decO,angleO,waveO,paO,timeLeftO,ra,dec,slitwidthO
integer,dimension(20) :: codeO
logical :: pastlog,exists
character(len=*),dimension(12),parameter :: month = (/ "jan", "feb", "mar", &
	"apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec" /)
character(len=6) :: white = "white"
!=====================================================
pastlog = ( iargc() > 0 )
if ( pastlog ) then
    call getarg(1,pastfile)
    datafile = TACTICAL//pastfile(:len_trim(pastfile))
    print*,datafile
    call plottargets()
else
    inquire(file=TACTICAL//"LOCK",exist=exists)   ! set LOCK
    if ( exists ) stop
    open(unit=99,file=TACTICAL//"LOCK",access="sequential",action="write", &
        status="new")
    datafile = TACTICAL//"telemetry1.txt"
    white = shade(1)
    do tick=0,900   ! each tick is "delay" seconds
        call soartelemetry()
        if ( modulo(tick,nint(15*60.0/delay)) == 0 ) call pasca()
        call plottargets()
        if ( tick > 0 .and. utnow >= sunrise ) exit
    end do
    write(unit=pastfile,fmt="(a,2i2.2,a)") "move "// &
        datafile(:len_trim(datafile))//" "//TACTICAL//month(mon),day, &
        yr-2000,"telem.txt"
    call system(pastfile)
    close(unit=99,status="scratch")       ! delete LOCK
end if

contains

subroutine soartelemetry()
character(len=180) :: cmd
integer :: hr,minute,deg,arcmin
real :: sec,arcsec,seeing,timUT,ra,dec
real :: f3,f4,f5,f8,f11,f14,f17,f18,f19,f21,f22,f23,f24,f25,f26,f27,f28,f29
integer :: i,i1,i2,i6,i7,i9,i10,i12,i13,i15,i16,statIO,ninst
character(len=10) :: c20,c30
character(len=3) :: stat
logical :: exists
character(len=*), parameter :: instrument = "AcqIFUGOOSPASOIOSIPHOVIS"
character(len=*), parameter :: f = "soartelescope.org/release/02about/eng_about/weather/telem.txt"
!=====================================================
cmd = WGETPATH//" -qr http://"//f
call system(cmd)
open(unit=1,file=f,access="sequential",action="read",status="old", &
	position="rewind",iostat=statIO)
read(unit=1,fmt="(a)") cmd
close(unit=1)
do i=11,len(trim(cmd))
	if ( cmd(i:i) == ":" ) cmd(i:i) = " "
end do
read(unit=cmd(11:),fmt=*,iostat=statIO) i1,i2,f3,f4,f5,i6,i7,f8,i9,i10,f11, &
        i12,i13,f14,i15,i16,f17,f18,f19,c20,f21,f22,f23,f24,f25,f26,f27, &
        f28,f29,c30
if ( statIO < 0 ) return
timUT = i1+(i2+f3/60.0)/60.0
hr = i6
minute = i7
sec = f8
ra = hr+(minute+sec/60.0)/60.0
deg = i9
arcmin = i10
arcsec = f11
if ( deg < 0 ) then
        dec = deg-(arcmin+arcsec/60.0)/60.0
else
        dec = deg+(arcmin+arcsec/60.0)/60.0
end if
seeing = f29
ninst = (index(instrument,c30(:3))-1)/3

datafile = TACTICAL//"telemetry1.txt"   ! need to do this to avoid corruption??
inquire(file=datafile,exist=exists)
if ( exists ) then
        stat = "old"
else
        stat = "new"
end if
print*,"#3 ",datafile
open(unit=1,file=datafile,access="sequential", &
        action="write",status=stat,position="append")
write(unit=1,fmt="(f7.4,3f8.3,f5.2,i3)") timUT,ra*15.0,dec,f5,seeing,ninst
close(unit=1)
end subroutine soartelemetry


subroutine pasca() 	! Process PASCA all-sky image into .bmp overlay for SN 
real, dimension(256,256) :: a
complex, dimension(256,256) :: data
integer, dimension(256,256) :: mask
character(len=120) :: cmd
real :: step,dmin,dmax, n2 = 18.0
integer :: nx,ny,nz,i,j
character(len=*), parameter :: f = &
"soartelescope.org/release/06observing/eng_observing/monitors/images/lastpix.fits"
!=====================================================
cmd = WGETPATH//" -qr http://"//f
call system(cmd)
call rfits(f,a,nx,ny,nz,dmin,dmax)
mask = 0
forall (i=1:nx, j=1:ny, (i-129)**2+(j-129)**2 <= 110**2 ) mask(i,j) = 1

call lineclean(a)
data = cmplx(mask*min(a,1000.0))
call four2(data,+1)			! remove residual stars
call filter(data,n2)			! low-pass filter
call four2(data,-1)
where ( mask > 0 )
	a = real(data)/(nx*ny)
elsewhere
	a = 0
end where
where ( mask > 0 ) a = a-minval(a, mask > 0)

step = maxval(a)/6.0		! set contour bins
a = step*ceiling(a/step)	! contour the clouds
a(1,1) = maxval(a)*1.5          ! this adjusts relative intensities
call filmod("delete")
call setfil("smlastpix.bmp")
call DSgraph("bmp")
call unit(0)
call IMage(a,parm="noxyz")
call disfin()
cmd = "imconvert smlastpix.bmp -trim "//TACTICAL//"current.bmp"
call system(cmd)                ! ImageMagick processing. trim & copy image
cmd = "mogrify -flip -geometry 950x950 -crop 950x950+70+60 "//TACTICAL//"current.bmp"
call system(cmd)                !ImageMagic scale/crop to fit SN display
end subroutine pasca

subroutine filter(a,rmax2)		! zero high frequency components (beyond rmax2)
complex, dimension(:,:), intent(inout) :: a
real, intent(in) :: rmax2
integer :: i,j,nx,ny
!=====================================================
nx = size(a,1)
ny = size(a,2)
do j=1,ny/2
	do i=1,nx/2
		if ( i**2+j**2 < rmax2 ) cycle
		a(i,j) = (0.0,0.0)
		a(nx+1-i,j) = (0.0,0.0)
		a(i,ny+1-j) = (0.0,0.0)
		a(nx+1-i,ny+1-j) = (0.0,0.0)
	end do
end do
end subroutine filter

subroutine lineclean(a) 	! simplified IRAF lineclean
real, dimension(:,:), intent(inout) :: a
integer, dimension(2), parameter :: ibc = (/ 3,3 /)
real, dimension(256) :: x,y,fit,diff
integer :: ix,iy,nx,i,mid
real :: ave,adev,sdev,var,skew,curt,splint
integer, parameter :: m = 15	! points to median (ODD!)
real, dimension(0:m-1) :: s
real, parameter :: sigmaclip = 1.5
!=====================================================
nx = size(a,1)
mid = nint(m/2.0+0.5)

do iy=1,size(a,2)
	i = 0
	do ix=1,nx-m,m
		i = i+1
		x(i) = ix+mid
		s = a(ix:ix+m-1,iy)
		call qsort(s)	! median
		y(i) = s(mid)
	end do
	call splred(x,y,diff,ibc,i,fit)
	fit(:m-1) = a(:m-1,iy)
	fit(nx-m:nx) = a(nx-m:nx,iy)
	do ix=m,nx-m
		fit(ix) = splint(real(ix),x,y,diff,i)
	end do
	diff = a(:,iy) - fit
	if ( any(diff /= 0.0) ) call moment(diff,ave,adev,sdev,var,skew,curt)
	where ( abs(diff) >= sdev*sigmaclip ) a(:,iy) = fit
end do
end subroutine lineclean


subroutine plottargets()
real :: alturaNOW
integer :: i
character(len=3) :: st
!=====================================================
call gettelem()
call getobjs()
call getview(epoch,yr,mon,day,utnow,object)
if ( pastlog ) then
    i = index(pastfile,"telem")
    read(unit=pastfile(i-7:i-1),fmt="(a,2i2)" ) st,day,yr
    yr = yr+2000
    do i=1,12
        if ( st == month(i) ) exit
    end do
    mon = i
    utnow = 9
else
    call gettimes(yr,mon,day,utnow)
end if
call addmoon(mon,day,utnow)
call precess(epoch,yr,mon,day)
call getnightime(sunset,tucrepv,tucrepm,sunrise)
call layout(utnow,sunset,tucrepv,tucrepm,sunrise,object)
call plotobs(utnow,sunset,tucrepv,tucrepm,sunrise,alturaNOW)
print*,object
if ( object > 0 ) call sliteff(tucrepv,slitWidthO(object),waveO(object),alturaNOW)
call disfin()
if ( WEB ) then
	write(unit=*,fmt="(a/)") "Content-type: image/gif"
	call system("/bin/cat ./dislin.gif")
!	call system("d:/unix/bin/cat --binary ./dislin.gif")
end if
end subroutine plottargets
	
subroutine sliteff(x,width,wave,alt)
real, intent(in) :: x,width,wave,alt
real :: see,wt
real, dimension(20), parameter :: mat = (/ 0.87, 0.77, 0.67, 0.59, 0.52, &
        0.95, 0.89, 0.82, 0.75, 0.68, 0.98, 0.95, 0.91, 0.85, 0.8, &
                0.99, 0.98, 0.96, 0.92, 0.88 /)
real, dimension(5), parameter :: w = (/ 0.45, 0.56, 0.85, 1.1, 1.5 /), &
        c = (/ 0.6, 0.8, 1.0, 1.2, 1.4 /)
character(len=30) :: cap
!===========================================================
if ( width <= 0.0 ) return
see = seeingAvg*(90.0-alt)**0.6*(550.0/wave)**1.2
print*,"object seeing, slitwidth, index :",alt,see,width,nint(see/0.2)-2
if ( see <= width .or. see <= 0.6 ) return
wt = -1
select case ( nint(see/0.2)-2 )
   case (1)
     if ( width <= 0.45 ) wt = 0.87
   case (2)
     if ( width <= 0.45 ) then
        wt = 0.77
     else if ( width <= 0.65 ) then
        wt = 0.89
     end if
   case (3)
      if ( width <= 0.45 ) then
         wt = 0.67
      else if ( width <= 0.65 ) then
          wt = 0.82
      end if
   case (4)
      if ( width <= 0.45 ) then
          wt = 0.59
      else if ( width <= 0.65 ) then
          wt = 0.75
      else if ( width <= 0.85 ) then
          wt = 0.85
      end if
   case (5)
      if ( width <= 0.45 ) then
          wt = 0.52
      else if ( width <= 0.65 ) then
          wt = 0.68
      else if ( width <= 0.85 ) then
          wt = 0.8
      else if ( width <= 1.1 ) then
          wt = 0.88
      end if
end select
if ( wt < 0.0 ) return
write(unit=cap,fmt="(a,i3,a)") "> ",nint((1-wt)*100.0)," % slit light loss"
if ( wt < 90.0 ) call color("yellow")
if ( wt < 75.0 ) call color("orange")
if ( wt < 0.6 ) call color("red")
call rlmess(cap,x*hr,80.0)
call color("white")
end subroutine sliteff

subroutine gettelem()
integer :: i,i6,statIO
real :: f1,f2,f3,f4,f5
integer, dimension(3) :: nSee
!=====================================================
open(unit=1,file=datafile,position="rewind", &
	action="read",status="old",iostat=statIO)
nC = 0
do i=1,2000
	read(unit=1,fmt=*,iostat=statIO) f1,f2,f3,f4,f5,i6
	if ( statIO < 0 .or. f1 <= 0.0 ) exit
	if ( f5 <= 0.0 ) cycle
	nC = nC+1
	utC(nC) = f1
	raC(nC) = f2       ! already degs
	decC(nC) = f3
        elevC(nC) = f4 
	seeing(nC) = f5
	inst(nC) = i6
end do
close(unit=1)
if ( nC >= 3 ) then
	where ( seeing(nC-2:nC) > 0.0 )
		nSee = 1
	elsewhere
		nSee = 0
	end where
	if ( sum(nSee) > 0 ) then
		seeingAvg = sum(seeing(nC-2:nC))/sum(nSee)
	else
		seeingAvg = 0.0
	end if
else
	seeingAvg = seeing(nC)
        nSee = 1
        if ( seeingAvg > 5.0 .or. seeingAvg < 0.01 ) seeingAvg = 0.0
end if
end subroutine gettelem

subroutine getobjs()
real :: arIN,decIN,timeleftIN,avgdec
integer :: i,ii,j1,j2,j3,j4,codeIN,paIN,waveIN,statIO,ic
character(len=200) :: line,nomIN
!=====================================================
open(unit=1,file=TACTICAL//"soar.txt",status="old",access="sequential", &
        position="rewind",action="read",iostat=statIO)
if ( statIO < 0 ) stop
read(unit=1,fmt="(a)",iostat=statIO) (line,i=1,21)	! skip 21|22 lines
if ( line(1:1) /= "1" ) read(unit=1,fmt="(a)",iostat=statIO) line
nO = 0
do
	read(unit=1,fmt="(a)",iostat=statIO) line
	if ( statIO < 0 .or. nO == 20 ) exit
	read(unit=line,fmt=*) j1,j2,j3,j4,arIN,decIN
	nO = nO+1
	raO(nO) = arIN/15.0	! to hours
	decO(nO) = decIN		! already degs
        avgdec = cos(0.5*(decIN+decC(nC))/57.29578) ! RA weight
        angleO(nO) = sqrt( (15.0*(arIN-raC(nC))*avgdec)**2 + (decIN-decC(nC))**2 )
        ic = index(line,"http")-1
        if ( ic < 0 ) then
           slitWidthO(nO) = 0.0
        else
           do i=ic,20,-1
               if ( line(i:i) == "\t" .or. line(i:i) == " " ) exit
           end do
           read(unit=line(i:ic),fmt=*) slitWidthO(nO)
        end if
	call getTarget(line,nomIN,codeIN,paIN,waveIN,timeleftIN)
	nomO(nO) = nomIN
	codeO(nO) = codeIN
	paO(nO) = paIN
	waveO(nO) = waveIN
	timeleftO(nO) = timeleftIN
	ii = max(index(line,"UNC"),index(line,"MSU"),index(line,"NOA"), &
		index(line,"BRA"),index(line,"SOA"))
	do i=max(1,ii),len(line)
		if ( line(i:i) == "\t" .or. line(i:i) == "_" .or. &
			line(i:i) == ":" ) line(i:i) = " "
	end do
	i = index(line,"http")-2
	if ( i <=0 ) i = 100
	ii = index(line,"f/")
	if ( ii <= 0 ) then
		msgO(nO) = " "
	else
		msgO(nO) = line(ii:i)
	end if
end do
close(unit=1)
end subroutine getobjs

subroutine getview(epoch,yr,mon,day,utnow,object)
integer,intent(out) :: epoch,yr,day,mon,object
real, intent(out) :: utnow
character(len=500) :: cgi
integer :: i1,i2,i,cmdlen
!=====================================================
if ( WEB ) then
	call getenv("QUERY_STRING",cgi)
!cgi = "http://localhost/cgi-bin/pngplot.cgi?wt=2000%2C%0D%0A0%2C0%2C0%2C%0D%0A0%2C%0D%0A1&but=PLOT"
	i1 = index(cgi,"wt=")+3
	i2 = index(cgi,"&but=")
	do i=i1,i2
		if ( cgi(i:i) == "+" ) then
			cgi(i:i) = " "
		else if ( cgi(i:i+2) == "%2C" ) then
			cgi(i:i+2) = " , "
		else if ( cgi(i:i+2) == "%0A" .or. cgi(i:i+2) == "%0D" ) then
			cgi(i:i+2) = "   "
		else if ( cgi(i:i+2) == "%2F" ) then
			cgi(i:i+2) = " / "
		end if
	end do
	cmdlen = 0
	do i=i1,i2
		if ( cgi(i:i) == "&" ) exit
		if ( cgi(i:i) /= " " ) then
			cmdlen = cmdlen+1
			cgi(cmdlen:cmdlen) = cgi(i:i)
		end if
	end do
	read(unit=cgi(:cmdlen),fmt=*) epoch,yr,day,mon,utnow !,object
else
	epoch = 0
	yr = 0
end if
object = minloc(angleO(:nO),dim=1) 
        ! confirm this w/ current instrument inst(nC) from telemetry
!print*,"#object,angle,instrument (should match that in use) ",object, &
!    object,codeO(object)
if ( angleO(object) > 5.0 ) object = 0
if ( epoch < 1950 .or. epoch > 2020 ) epoch = 2000
end subroutine getview

subroutine getTarget(L,N,codeIN,paIN,waveIN,timeleft)
character(len=140),intent(in) :: L
character(len=30),intent(out) :: N
integer,intent(out) :: codeIN,paIN,waveIN
real,intent(out) :: timeleft
real :: time
integer :: i,ip,j
! 1 = white, 2 = red, 3 = green, 4 = blue, 5 = cyan, 6 = pink, 7 = yellow,
! 8 = orange, 9 = brown
character(len=*),parameter :: inst = "...PHXSPN...GSPVIMIFUOSIWFSVIS"
! code = 2 (PHX), 3 (SPN), 5(GSP), 6(VIM), 7(IFU), 8(OSI), 9(CAL)
!these set the order of color codes in SOARstatus1 & control1
!===================================================
ip = index(L,"IFU")+index(L,"GSP")+index(L,"SPN")+index(L,"VIM")+ &
	index(L,"OSI")+index(L,"PHX")+index(L,"CAL")+index(L,"WFS")+ &
	index(L,"VIS")
if ( ip <= 0 ) then
!	      print*,"Invalid object"
	return
else
	codeIN = max(0,(index(inst,L(ip:ip+2))-1)/3+1)
end if
! eat up blanks after instrument name
do i=ip+3,len(L)
	if ( L(i:i) /= " " .and. L(i:i) /= "\t" ) exit
end do
N = " "
j = 0
do ip=i,len(l)	! eat up target name which ends w/ a |
	if ( L(ip:ip) == "|" .or. L(ip:ip) == "\t" ) exit
	j = j+1
	if ( L(ip:ip) == "_" ) then
		N(j:j) = " "
	else
		N(j:j) = L(ip:ip)
	end if
end do
ip = ip+1
if ( index(L(ip:),"Any")+index(L(ip:),"f/") > 0 ) then
	read(unit=L(ip:),fmt=*) time,timeleft,waveIN
	paIN = 0.0
else
	read(unit=L(ip:),fmt=*) time,timeleft,waveIN,paIN
end if
if ( waveIN <= 300.0 .or. waveIN >= 8000.0 ) waveIN = 550.0
!print*,'pa=',paIN,' degs. Wavelength=',waveIN,' nm. Timeleft=',timeleft
end subroutine getTarget

subroutine gettimes(yr,mon,day,utnow)
integer,intent(inout) :: yr,mon,day
real, intent(inout) :: utnow
character(len=32) :: dd,timer,timezone
integer :: hh,mm,ss
!=====================================================
call date_and_time(dd,timer,timezone)
! some sanity checks
if ( yr < 2006 .or. yr > 2020 .or. mon < 1 .or. mon > 12 .or. day < 1 &
	.or. day > 31 .or. .not.WEB ) read(unit=dd,fmt="(i4,2i2)") yr,mon,day
if ( .not.WEB .or. utnow <= 0.0 .or. utnow > 24.0 ) then
	read(unit=timer,fmt="(3i2)") hh,mm,ss
	utnow = hh+(mm+ss/60.0)/60.0 -19
	if ( utnow >= 24.0 ) then
		utnow = utnow - 24.0
		day = day+1
	else if ( utnow < 0.0 ) then
		utnow = utnow + 24.0
	end if
        if ( utnow > sunset ) day = day+1
end if
end subroutine gettimes

subroutine addmoon(mon,day,utnow)
integer,intent(inout) :: mon,day
real,intent(in) :: utnow
!=====================================================
nO = nO+1
codeO(nO) = 1
nomO(nO) = "Moon (@ UT now)"
call moonp1(long,glat,yr,mon,day,utnow,ratop,dectop)
raO(nO) = ratop
decO(nO) = dectop
end subroutine addmoon

subroutine precess(epoch,yr,mon,day)
integer,intent(in) :: epoch,yr,mon,day
integer :: month,i
real :: tff,ti,tf,gio,zeta,teta
real,dimension(3,3) :: m
real, dimension(3) :: x0,x
!=====================================================
if (mon <= 2) then
	month = aint((mon-1)*63.0/2.0)
else
	month = aint((mon+1)*30.6)-63
end if
month = month+day
tfF = yr+month/365.25
ti = (epoch-2000.0)/100.0
tf = (tfF-2000.0-100.0*ti)/100.0
! elementos precesionales en grados
gio = tf*(2306.2181+1.39656*ti-0.000139*ti**2+(0.30188-0.000344*ti)*tf &
	+0.017998*tf**2)/3600.0
zeta = gio+tf**2*(0.7928+0.00041*ti+0.000205*tf)/3600.0
teta = tf*(2004.3109-0.8533*ti-0.000217*ti**2-(0.42665+0.000217*ti)*tf &
	-0.041833*tf**2)/3600.0
! convert to rads
gio = gio/r2d
zeta = zeta/r2d
teta = teta/r2d
! define rotation matrix
M(1,1)=-(sin(gio)*sin(zeta))+cos(gio)*cos(teta)*cos(zeta)
M(1,2)=-(cos(gio)*sin(zeta))-(sin(gio)*cos(zeta)*cos(teta))
M(1,3)=-sin(teta)*cos(zeta)
M(2,1)=sin(gio)*cos(zeta)+cos(gio)*cos(teta)*sin(zeta)
M(2,2)=cos(gio)*cos(zeta)-sin(gio)*cos(teta)*sin(zeta)
M(2,3)=-sin(teta)*sin(zeta)
M(3,1)=cos(gio)*sin(teta)
M(3,2)=-sin(gio)*sin(teta)
M(3,3)=cos(teta)

do I=1,nO
	X0(1) = cos(decO(I)/r2d)*cos(raO(I)*15.0/r2d)
	X0(2) = cos(decO(I)/r2d)*sin(raO(I)*15.0/r2d)
	X0(3) = sin(decO(I)/r2d)
	x = matmul(m,x0)
	ra(I) = r2d*ATAN2(X(2),X(1))/15.0
	if ( ra(I) < 0.0 ) ra(I) = ra(I)+24.0
	dec(I) = r2d*Asin(X(3))
	if ( I == nO ) then	! moon
		ra(I) = ratop
		dec(I) = dectop
	end if
end do
end subroutine precess

subroutine getnightime(sunset,tucrepv,tucrepm,sunrise)
real,intent(out) :: sunset,sunrise,tucrepv,tucrepm
real :: rearth,anghoa,diszen,decsun,arsun,hocaso,tucrepm1, &
	tuocaso1,tuocaso2,tucrepv1,tucrepv2,tuorto1,horto,tuorto2,tucrepm2
! DISZEN: refraccion 34'
!         semidiametro del Sol: 16'
!         depresion del horizonte: usamos elipsoide WGS-84
real,parameter :: semia = 6378.137, semib = 6356.752
!=====================================================
! Calculate Julian date (FJ) and TS at 0h TU Greenwich(ts0)
call fj_ts0(yr,mon,day,0.0)
!                                     POSICION DEL SOL, ORTO Y OCASO
! Iteramos para calcular la posicion del Sol en el instante del ocaso y
! del orto. Tomamos como valor inicial la posicion a 0h T.U. del dia
! de la observacion.
call POSSOL(0,arsun,decsun)
REARTH=1.0/SQRT( cos(LAT/r2d)*cos(LAT/r2d)/SEMIA**2 + &
	sin(LAT/r2d)*sin(LAT/r2d)/SEMIB**2 )
ANGHOA=r2d*acos(REARTH/(REARTH+ALTOBS/1000.0))
DISZEN=90.0+34.0/60.0+16.0/60.0+ANGHOA
call ORCASO(DECSUN,DISZEN,horto,hocaso)
! ocaso del Sol (calculamos el ocaso que tiene lugar a partir de las 0h
! de T.U. del dia de la observacion).
call tiEMPOU(ARSUN,HOCASO,LONG,ts0,TUOCASO1)
call test(0,tuocaso1,diszen,arsun,decsun,horto,hocaso,hocaso,tuocaso2)
! crepusculo astronomico vespertino (utilizamos como valor inicial en
! la iteracion el instante del ocaso)
tucrepv1 = tuocaso2
call test(0,tucrepv1,108+anghoa,arsun,decsun,horto,hocaso,hocaso,tucrepv2)
! orto del Sol (calculamos el orto que tiene lugar a partir de las 0h
! de T.U. del dia de la observacion)
tuorto1 = tuocaso2
call test(0,tuorto1,diszen,arsun,decsun,horto,hocaso,horto,tuorto2)
! crepusculo astronomico matutino (utilizamos como valor inicial en
! la iteracion el instante del orto)
tucrepm1 = tuorto2
call test(0,tucrepm1,108+anghoa,arsun,decsun,horto,hocaso,horto,tucrepm2)
! si el orto calculado es anterior al ocaso quiere decir que tenemos
! que recalcular el orto que ocurre para el dia siguiente al dia de la observacion.
if (TUORTO2 < TUOCASO2) then
	TUORTO1=TUORTO2
	call test(1,tuorto1,diszen,arsun,decsun,horto,hocaso,horto,tuorto2)
! crepusculo astronomico matutino (utilizamos como valor inicial en
! la iteracion el instante del crepusculo para el dia anterior)
	TUCREPM1=TUCREPM2
	call test(1,tucrepm1,108+anghoa,arsun,decsun,horto,hocaso,horto,tucrepm2)
end if

! limites temporales en los cuales trabajaremos
!print*,"#1 sunset,tucrepv,tucrepm.sunrise:",tuocaso2,tucrepv2,tucrepm2,tuorto2
sunset = TUOCASO2 +15
if ( sunset > 24.0 ) sunset = sunset-24.0
TUCREPV = TUCREPV2 +15
if ( tucrepv > 24.0 ) tucrepv = tucrepv-24.0
TUCREPM = TUCREPM2 +15
if ( tucrepm > 24.0 ) tucrepm = tucrepm-24.0
sunrise = TUORTO2 +15
if ( sunrise > 24.0 ) sunrise = sunrise-24.0
!if (TUCREPV < sunrise) TUCREPV=TUCREPV+24
!print*,"#2 sunset,tucrepv,tucrepm.sunrise:",sunset,tucrepv,tucrepm,sunrise
! recalculamos la FJ y ts0 para el dia de la observacion a 0h TU
call fj_ts0(yr,mon,day,0.0)
end subroutine getnightime

subroutine POSSOL(IAP,arsun,decsun)
real,intent(out) :: arsun,decsun
!	   Posicion del Sol, J. Meeus
!	   Astronomical Formulae for Calculators
!	   Cap. 18, pag. 79 y siguientes
!          Subrutinas escritas por M. Cornide
! IAP = 1 ---> coordenadas aparentes
! IAP = 0 ---> coordenadas no aparentes
integer,intent(in) :: IAP
real(kind=double) :: ar,ALAMBDA,PAR,DEC,zero = 0
!=====================================================
call PSOL(FJ,IAP,ALAMBDA,PAR)
call ECUAT(FJ,ALAMBDA,zero,AR,DEC)
ARSUN = AR*r2d/15.0
DECSUN = DEC*r2d
end subroutine possol

subroutine PSOL(FJ,IAP,ALON,PAR)
real(kind=double),intent(out) :: alon,par
!	      AL : Long. media
!	      AM : Anomalia media
!	      EXC: Excentricidad
!	      EC : Ecuacion del centro
!	      V  : Anomalia verdadera
!	      R  : Radiovector
real(kind=double) :: t,al,am,exc,c1,c2,ec,v,r,a,b,c,d,e,h,om,rr
real(kind=double),intent(in) :: FJ
real(kind=double),parameter :: DPI=2*PI, full = 360
integer,intent(in) :: iap
!=====================================================
T = (FJ-2415020.0)/36525.0
AL = 279.69668+36000.76892*T+0.0003025*T*T
AM = 358.47583+35999.04975*T-0.000150*T*T+0.0000033*T*T*T
AL = modulo(AL,full)/r2d
AM = modulo(AM,full)/r2d
EXC = 0.016751014-0.0000418*T-0.000000126*T*T
 C1 = 1.919460-0.004789*T-0.000014*T*T
 C2 = 0.020094-0.000100*T
EC = (C1*sin(AM)+C2*sin(2*AM)+0.000293*sin(3*AM))/r2d

ALON = AL+EC
V = AM+EC
R = 1.0000002*(1.0-EXC*EXC)/(1.0+EXC*cos(V))
!    Perturbaciones
A = (153.23+22518.7541*T)/r2d
B = (216.57+45037.5082*T)/r2d
C = (312.69+32964.3577*T)/r2d
D = (350.74+445267.1142*T-0.00144*T*T)/r2d
E = (231.19+20.20*T)/r2d
H = (353.40+65928.7155*T)/r2d

A = modulo(A,DPI)
B = modulo(B,DPI)
C = modulo(C,DPI)
D = modulo(D,DPI)
H = modulo(H,DPI)
E = modulo(E,DPI)

ALON = ALON+(0.00134*cos(A)+0.00154*cos(B)+0.00200*cos(C)+ &
	0.00179*sin(D)+0.00178*sin(E))/r2d
R = R+0.00000543*sin(A)+0.00001575*sin(B)+0.00001627*sin(C)+ & 
	0.00003076*cos(D)+0.00000927*sin(H)

if ( IAP == 1 ) then
	OM = (259.18-1934.142*T)/r2d
	OM = modulo(OM,DPI)
	if ( OM < 0.0) OM = OM+DPI
	ALON = ALON-(0.00569+0.00479*sin(OM))/r2d
end if
RR = R*1.49597870e11/6.37814e6
PAR = asin(1.0/RR)*r2d
end subroutine psol

subroutine EPSIL(FJ,EPSI)
real(kind=double),intent(in) :: FJ
real(kind=double),intent(out) :: epsi
real(kind=double) :: u,a1,a2
!=====================================================
U = (FJ-2451545.0)/3652500.0
!  ----- Calculo de la oblicuidad de la ecliptica
A1 = 2.18-3375.70*U+0.36*U**2
A2 = 3.51+125666.39*U+0.10*U**2
EPSI = 0.4090928-0.0226938*U-75.0e-7*U**2+96926.0e-7*U**3 &
    -2491.0e-7*U**4-12104.0e-7*U**5+(446*cos(A1)+28*cos(A2))*1.0e-7
end subroutine epsil

subroutine ECUAT(FJ,AL,B,A,D)
real(kind=double),intent(out) :: a,d
real(kind=double),intent(in) :: FJ
real(kind=double),intent(in) :: al,b
real(kind=double) :: x0,y0,z0,x,y,z,eps
!	    Calcula  Ascension Recta y Declinacion
real(kind=double), parameter :: DPI = 2*PI
!=====================================================
call EPSIL(FJ,EPS)
X0 = cos(B)*cos(AL)
Y0 = cos(B)*sin(AL)
Z0 = sin(B)
X = X0
Y = Y0*cos(EPS)-Z0*sin(EPS)
Z = Y0*sin(EPS)+Z0*cos(EPS)
A = atan2(Y,X)
if ( A < 0.0 ) A = A+DPI
D = asin(Z)
end subroutine ecuat

subroutine fj_ts0(yr,mon,day,hour)
integer,intent(in) :: yr,mon,day
real,intent(in) :: hour
! **********************************************************************
!                                                      SUBROUtiNE FJ_ts0
!                                                      *****************
! Calcula la fecha juliana y el Tiempo Sidereo en Greenwich a 0h UT a
! partir de los datos de una fecha concreta. En la fecha juliana se
! tiene en cuenta la fraccion del dia correspondiente a la hora
! traspasada a traves de la variable global HORA.
! Los calculos se realizan en doble precision para no perder informacion
! por redondeo.
integer :: FECHA
real(kind=double) :: M,A,B,Y,DT,Dts0,DANO,DMES,DDIA
real(kind=double) :: v = 8.64e-4
!=====================================================
DANO = yr
DMES = mon
DDIA = day
FECHA = int(yr)*10000+int(mon)*100+int(day)
if ( FECHA >=  15821015 ) then
	A = int(DANO/100.0)
	B = 2.0-A+int(A/4.0)
else
	B = 0.0
end if
if ( yr >=  0.0 ) then
	A = 0.0
else
	A = -0.75
end if
if ( mon >=  2.0 ) then
	Y = DANO
	M = DMES
else
	Y = DANO-1.0
	M = DMES+12.0
end if
FJ = int(365.25*Y+A)+int(30.6001*(M+1.0))+DDIA+B+1720994.5
DT = (FJ-2451545.0)/36525.0
Dts0 = 6.0*3600.0+41.0*60.0+50.54851+8640184.812866*DT &
	+0.093104*DT*DT-6.2e-6*DT*DT*DT
Dts0 = modulo(Dts0,v)
Dts0 = Dts0/8.64e4
Dts0 = Dts0*24.0
ts0 = Dts0
if ( ts0 < 0.0 ) then
	ts0 = ts0+24.0
else if ( ts0 > 24.0 ) then
	ts0 = ts0-24.0
end if
FJ = FJ + hour/24.0
end subroutine fj_ts0

subroutine orcaso(dec,diszen,horto,hocaso)
real,intent(in) :: dec,diszen
real,intent(out) :: horto,hocaso
! **********************************************************************
!                                                    SUBROUtiNE ORCASO
!                                                      *****************
! Calculamos para un objeto concreto, su visibilidad, ortos, ocasos, 
! culminacion y condiciones en las cuales obtiene una distancia
! cenital igual a DISZEN.
REAL :: LATR,DECR,COSENOH,COSENOA,ALTR

REAL :: AORTO,AOCASO,ALTCUL
CHARACTER(len=1) :: ESTADO,DONCUL      
!=====================================================
! Determinamos visibilidad de los objetos y determinamos su ESTADO
! segun el siguiente criterio
!     ESTADO = I ==> Objeto INVISIBLE
!     ESTADO = C ==> Objeto CIRCUMPOLAR
!     ESTADO = O ==> Objeto con ORTOS y OCASOS
if ( LAT >= 0 ) then
	if ( DEC >= 90.0-LAT ) then
          ESTADO="C"
	else if ( DEC <= LAT-90.0 ) then
          ESTADO="I"
	else
          ESTADO="O"
	end if
else
	if ( DEC <= -(90.0+LAT) ) then
          ESTADO="C"
	else if ( DEC >= 90.0+LAT ) then
          ESTADO="I"
	else
          ESTADO="O"
	end if
end if
LATR=LAT/r2d
DECR=DEC/r2d
ALTR=(90-DISZEN)/r2d
if ( ESTADO == "O" ) then
	COSENOH=(sin(ALTR)-sin(LATR)*sin(DECR))/(cos(LATR)*cos(DECR))
	if ( cosenoh < -1.0) cosenoh=-1
	if ( cosenoh > +1.0) cosenoh=+1
	COSENOA=(sin(ALTR)*sin(LATR)-sin(DECR))/(cos(ALTR)*cos(LATR))
	if ( cosenoa < -1.0) cosenoa=-1
	if ( cosenoa > +1.0 ) cosenoa=+1
	HOCASO=acos(COSENOH)
	HOCASO=HOCASO*r2d
	AOCASO=acos(COSENOA)
	AOCASO=AOCASO*r2d
	HORTO=360-HOCASO
	AORTO=360-AOCASO
end if
HORTO=HORTO/15.0
HOCASO=HOCASO/15.0
! Determinamos la altura en las culminaciones superiores
if ( ESTADO /= "I" ) then
     if ( LATR >= 0 ) then
          if ( DEC <= LAT ) then
            ALTCUL=90-(LATR-DECR)*r2d
            DONCUL="S"
          else
            ALTCUL=90-(DECR-LATR)*r2d
            DONCUL="N"
          end if
     else
          if ( DEC <= LAT ) then
            ALTCUL=90-ABS(DECR-LATR)*r2d
            DONCUL="S"
          else
            ALTCUL=90-ABS(LATR-DECR)*r2d
            DONCUL="N"
          end if
     end if
end if
end subroutine orcaso

subroutine tiempou(A,H,L,ts0,TU)
! Calculo del Tiempo Universal.
! Datos: A.R., h, Longitud, TSG0
real,intent(in) :: A,H,L,ts0
real,intent(out) :: TU
real :: TSL,TSG,DTS
!=====================================================
TSL = A+H
if ( TSL > 24.0 ) TSL = TSL-24.0
TSG = TSL-L
if ( TSG > 24.0 ) TSG = TSG-24.0
if ( TSG < 0.0 ) TSG = TSG+24.0
DTS = TSG-ts0
if ( DTS < 0.0 ) DTS = DTS+24.0
TU = 0.997270*DTS
if ( TU < 0.0 ) TU = TU+24.0
end subroutine tiempou

subroutine layout(utnow,sunset,tucrepv,tucrepm,sunrise,object)
real, intent(inout) :: utnow
integer,intent(in) :: object
real,intent(inout) :: sunset,tucrepv,tucrepm,sunrise
real, dimension(2) :: x,y
real,dimension(8),parameter :: val = (/1.0,1.1,1.2,1.3,1.4,1.5,2.0,3.0 /)
integer :: i,hours,mins
character(len=80) :: cap
real :: lo,tmp,first
!=====================================================
if ( sunset > sunrise ) sunset = sunset-24
if ( tucrepv > sunrise ) tucrepv = tucrepv-24
!print*,"utnow,sunset,tucrepv,tucrepm,sunrise:",utnow,sunset,tucrepv,tucrepm,sunrise
if ( WEB .or. pastlog ) then
        write(unit=cap,fmt="(a,2i2.2,a)") "c:/Apache/htdocs/"//month(mon), &
                day,yr-2000,".png"
        call setfil(cap)
	call metafl("png")
	if ( WEB) call filmod("delete")
else
	call metafl("xwin")
end if
call setpag("DA4P")
call disini()
call errmod("PROTOCOL","OFF")
call winmod("delay")
call winopt(delay,"delay")         ! delay til next update
call unit(0)
call sclmod("full")
call pagfll(4)
call height(30)
call axspos(250,2730)
call axslen(1700,2200)
if ( pastlog ) utnow = tucrepm
call color(white)
call messag("UT",125,2770)
write(unit=cap,fmt="(a,i3.2)") month(mon),day
call messag(cap,10,2690)
call frame(1)
call setgrf("name","none","none","none")
call labdig(-1,"x")
call labdig(-1,"y")
call labels("hours","x")
!print*,"#4,sunset,sunrise",sunset,sunrise
if ( sunset < 0.0 ) then
        first = 0.0
else
        first = int(sunset+1.0)
end if
call graf(sunset*hr,sunrise*hr,first*hr,hr,10.0,90.0,10.0,10.0)
do i=20,80,10
	call rlnumb(real(i),-1,(sunrise-0.5)*hr,real(i))
end do
call dot()
call color("blue")
x(1) = sunset*hr
x(2) = sunrise*hr
do i=10,80,10
	y = i
	call curve(x,y,2)
end do
y(1) = 10.0
y(2) = 90.0
do i=nint(sunset)+1,nint(sunrise)+12
	x = i*hr
	call curve(x,y,2)
end do
call dash()
call color(white)
x = sunrise*hr
call rlmess("sunrise",x(1)-0.5*hr,92.0)
if ( sunset > 20.0 ) then
	x = (sunset-24)*hr
else
	x = sunset*hr
end if
call rlmess("sunset",x(1)-0.5*hr,92.0)
if ( tucrepv > 20.0 ) then
	x = (tucrepv-24)*hr
else
	x = tucrepv*hr
end if
lo = x(1)
call rlmess("end twilight",x(1)-0.5*hr,95.0)
call curve(x,y,2)
x = tucrepm*hr
call rlmess("start dawn",x(1)-0.5*hr,95.0)
call curve(x,y,2)
call height(40)
if ( .not.pastlog ) then
     if ( utnow > 20.0 ) then
	    x = (utnow-24)*hr	! dislin silliness
    else
	    x = utnow*hr
    end if
    call color("red")
    call curve(x,y,2)
    call rlmess("NOW",x(1),89.0)
    call color("yellow")
    if ( utnow > lo/hr .and. utnow < tucrepm ) then
	    hours = tucrepm-utnow
	    mins = nint( (tucrepm-utnow-hours)*60.0 )
	    write(unit=cap,fmt="(i2,a,i2.2,a)") hours,":",mins," hrs until twilight"
    else
	    write(unit=cap,fmt="(f5.1,a)") tucrepm-lo/hr," hrs total dark"
    end if
    tmp = tucrepv
    if ( tucrepv > 24.0 ) tmp = tucrepv-24.0
    call rlmess(cap,tmp*hr,86.0)
    if ( utnow > sunset .and. utnow < sunrise ) then
	    hours = sunrise-utnow
	    mins = nint( (sunrise-utnow-hours)*60.0 )
	    write(unit=cap,fmt="(i2,a,i2.2,a)") hours,":",mins," hrs until sunrise"
	    call rlmess(cap,tmp*hr,83.0)
    end if
    call color(white)
    call solid()
    if ( seeingAvg > 0.1 ) then
        do i=1,8
                tmp = seeingAvg*val(i)**0.6
                if ( object > 0 ) tmp = tmp*(550.0/waveO(object))**1.2
	        call makecap(90-ACOS(1.0/val(i))*r2d,sunset*hr,sunrise*hr, &
		        tmp,val(i))
        end do
    end if
end if
call solid()
call height(26)
do i=1,nO-1
	if ( i == object ) then
		call color("red")
	else
                call color(shade(codeO(i)))
	end if
	write(unit=cap,fmt="(i2,t4,a,i5,f6.2,i5,a)") i,nomO(i), &
		nint(waveO(i)),timeleftO(i),nint(paO(i))," "//msgO(i)
	call messag(cap,20,45*i)
end do
call color(white)
call plotseeing(sunset+1.0)
end subroutine layout

subroutine makecap(y,tuini,tufin,seeing,val)
real,intent(in) :: y,seeing,tuini,tufin,val
!=====================================================
call height(26)
call rlnumb(val,1,tufin,y)
call height(50)
call rlnumb(seeing,1,tuini-1.12*hr,y)
end subroutine makecap

subroutine plotobs(utnow,sunset,tucrepv,tucrepm,sunrise,alturaNOW)
real,intent(in) :: utnow,sunset,tucrepv,tucrepm,sunrise
real, intent(out) :: alturaNOW
integer :: i,k
logical :: moontrack
real :: hour,tsg,tsl,hor
real,dimension(NP) :: tu,acimut,altura
!=====================================================
call height(20)
do i=1,nO
	moontrack = ( i == nO )
	do k=1,NP
		tu(k) = sunset+(sunrise-sunset)*(k-1.0)/(NP-1)
		if ( tu(k) >= 24.0 ) tu(k) = tu(k)-24.0
		TSG=ts0+TU(K)*1.002737909
! TSL para el instante TU
		TSL=TSG+LONG
		if ( TSL >= 24.0 ) TSL = TSL-24.0
! angulo horario
		if ( moontrack ) then
			hour = tu(k)
			call moonp1(long,glat,yr,mon,day,hour,ratop,dectop)
			HOR=TSL-ratop+0.5
		else
			HOR=TSL-ra(i)
			dectop = dec(i)
		end if
                HOR = HOR+8
		if ( hor >= 24.0 ) then
                        hor = hor-24.0
                        day = day+1
                end if
! pasamos a coordenadas horizontales
		call CAMBCOOR(HOR,dectop,ACIMUT(K),ALTURA(K))
	end do
    if ( pastlog .and. codeO(i) == 1 ) then
        call color(white)
    else
	    call color(shade(codeO(i)))
    end if
    write(unit=cap,fmt="(i2)") i
    call dibuobj(i,cap,hor,utnow,moontrack,tu,acimut,altura)
    if ( i == 1 ) alturaNOW = altura(1)
end do
call height(36)
end subroutine plotobs

subroutine CAMBCOOR(HOR,DEC,ACI,ALT)
! Transformaciones de coordenadas:
! (angulo horario,declinacion) ----> (acimut,altura)
real,intent(in) :: HOR,DEC
real,intent(out) :: ACI,ALT
real :: HORR,DECR,LATR,SENOAC,cosEAC
!=====================================================
HORR = HOR*15.0/r2d
DECR = DEC/r2d
LATR = LAT/r2d
SENOAC = cos(DECR)*sin(HORR)
cOSEAC = -(sin(DECR)*cos(LATR))+cos(DECR)*sin(LATR)*cos(HORR)
ACI = r2d*ATAN2(SENOAC,cosEAC)
if ( ACI < 0.0 ) ACI = ACI+360.0
ALT = r2d*Asin(sin(LATR)*sin(DECR)+cos(LATR)*cos(DECR)*cos(HORR))
end subroutine cambcoor


subroutine moonp1(long,lat,yr,mon,day,hour,ratop,dectop)
integer,intent(in) :: yr,mon,day
real,intent(out) :: ratop,dectop
real,intent(in) :: long,lat,hour
real(kind=double), parameter :: zero = 0, one = 1
!=====================================================
call moonp(one*15*long,one*lat,yr,mon,one*day,one*hour,zero, &
	ratop,dectop)
end subroutine moonp1

subroutine dibuobj(target,cap,HOR,utnow,moontrack,tu,acimut,altura)
integer,intent(in) :: target
character(len=2),intent(in) :: cap
real,intent(in) :: utnow,hor
logical,intent(in) :: moontrack
real,intent(in),dimension(NP) :: tu,acimut,altura
real :: parangle,utdiff,ut,x,y,x01,y01
integer :: j,show
logical :: graf, first
!=====================================================
first = .true.
show = 0
do j=1,NP
	graf = ( ALTURA(j) >= 10.0 )
	if ( .not.graf ) cycle
	Y=ALTURA(j)
	ut = tu(j)
	if ( ut >= 24.0 ) ut = ut-24.0
	X=ut*hr
	if ( ut < 0.0 ) ut = ut+24.0
	utdiff = ut-utnow
	call linwid(1)
	if ( .not.moontrack .and. utdiff > 0.0 .and. codeO(target) == 5 ) &
		call SetPA(target,utdiff,altura(j),object,1,parangle)
	if ( first ) then
		first = .false.
	else
		call rline(x01,y01,x,y)
		if ( .not.moontrack .and. modulo(j,20) == 0 ) call rlmess(cap,x,y)
	end if
	X01=X
	Y01=Y
! report parallactic angle if within 2 hrs of UT Now & if spectrometer
	if ( moontrack .or. utdiff < 0.0 .or. utdiff > 2.0 .or. codeO(target) /= 5 ) cycle
	call height(20)
	call color(shade(codeO(target)))
	call SetPA(target,utdiff,altura(j),object,2,parangle)
	if ( modulo(show,4) == 0 ) then ! reduce visual clutter
		show = 1			! show every 3rd
		call linwid(1)
		call rlnumb(parangle,-1,x,y-2)
	else
		show = show+1
	end if
end do
end subroutine dibuobj

subroutine SetPA(target,utdiff,altura,object,id,parangle)
real,intent(in) :: utdiff,altura
integer,intent(in) :: object,id,target
real,intent(out) :: parangle
real :: aa
!================================================== present & future
aa = sin(utdiff*15/r2d)*sin((90-lat)/r2d)/sin((90-altura)/r2d)
parangle = modulo(r2d*asin(max(-1.0,min(aa,1.0))),180.0)
if ( parangle < 0.0 ) parangle = parangle+180
! PA must agree w/ parallactic to within +/- 15 deg, for either slit orientation
if ( abs(paO(object)-parangle) <= 15.0 .or. abs(180+paO(object)-parangle) <= 15.0 ) then
	if ( id == 1 ) then
		call linwid(18) ! thick line
	else 
		call height(25) ! large PA caption
	end if
	if ( timeLeftO(object) >= utdiff ) then
		call color("red")   ! valid interval is red
	else			    ! full exposure is thick line
		call color(shade(codeO(target)))
	end if
end if
end subroutine SetPA

subroutine plotseeing(start)
real,intent(in) :: start
integer :: i,nxposn
!=====================================================
call height(28)
call labdig(1,"y")
call yaxis(0.0,1.5,0.0,0.2,1100," ",1,nxposn(start*hr),2730)
call hsymbl(8)
do i=1,nC
	if ( utC(i) >= start .and. utnow >= utC(i) &
		.and. seeing(i) > 0.1 ) then
	        call rlsymb(21,utC(i)*hr,10.0+seeing(i)*40.0/1.5)
        end if
end do
end subroutine plotseeing

subroutine test(ind,a,c,arsun,decsun,horto,hocaso,d,b)
integer,intent(in) :: ind
real,intent(inout) :: a,b,horto,hocaso,arsun,decsun,d
real,intent(in) :: c
real :: hora
integer :: i
!=====================================================
do i=1,4
	if ( ind > 0 ) then
		hora = a+24
	else
		hora = a
	end if
	call fj_ts0(yr,mon,day,hora)
	if ( ind > 0 ) then
		ts0 = ts0 + 24.0*2.737909e-3
		if ( ts0 >= 24.0 ) ts0 = ts0-24
	end if 
	call POSSOL(0,arsun,decsun)
	call ORCASO(DECSUN,c,horto,hocaso)
	call tiEMPOU(ARSUN,d,LONG,ts0,b)
	if ( b < 0.0 ) b = b+24
	if ( abs(a-b) <= 0.0003 ) exit
		a=b
end do
if ( b > 24.0 ) b = b-24
end subroutine test

function FNday(y,m,d,h) result(r)
!   FNday only works between 1901 to 2099 - see Meeus chapter 7
real(kind=double) ,intent(in) :: d,h
integer,intent(in) :: y,m
real(kind=double) :: r
!=====================================================
r = 367*y - 7*(y + (m+9)/12)/4 + (275*m)/9 + d - 730531.5 + h/24
end function Fnday

function FNrange(x) result(r)
!   the function below returns the true integer part, even for -ve numbers
real(kind=double),intent(in) :: x
real(kind=double) :: r
real(kind=double),parameter :: twoPi = 2*PI
!=====================================================
r = x - int(x / twopi) * twopi
end function FNrange


subroutine moonp(glong,glat,y,m,day,hr,mins,ratop,dectop)
!   Moon positions to a quarter of a degree on the sky
!   From page D76 of Astronomical Almanac
!   max RA error 97 (seconds of time), rms error 22 secs
!   max dec error 811 arcseconds, rms error 224 arcseconds
! compared to Integrated Computer Ephemeris for 18.6 years either side of J2000
!       program moontest
!	call moon(-1.91667d0,52.5d0,1998,8,9d0,11d0,56d0,ra,dec)
!        print*,'ra,dec = ',ra,dec
!        end program moontest
real(kind=double),intent(in) :: glong,glat,day,hr,mins
real(kind=single),intent(out) :: ratop,dectop
real(kind=double) :: long,lat,ra,dec,d,h,twoPI,degs,rads,t,l,bm, &
	gp,sdia,rm,xg,yg,zg,ecl,xe,ye,ze,lst,xtop,ytop,ztop,rtop
integer,intent(in) :: y,m
!=====================================================
!        print*,'input: ',glong,glat,y,m,day,hr,mins
twopi = 2*PI
degs = r2d
rads = 1.0/r2d
h = hr + mins / 60.0
d = FNday(y, m, day, h)
t = d / 36525.0
!
!   calculate the geocentric longitude
l = FNrange(rads * (218.32 + 481267.883 * t))
l = l + rads*6.29 * sin(FNrange((134.9 + 477198.85 * t)*rads))
l = l - rads*1.27 * sin(FNrange((259.2 - 413335.38 * t)*rads))
l = l + rads*0.66 * sin(FNrange((235.7 + 890534.23 * t)*rads))
l = l + rads*0.21 * sin(FNrange((269.9 + 954397.7 * t)*rads))
l = l - rads*0.19 * sin(FNrange((357.5 + 35999.05 * t)*rads))
l = l - rads*0.11 * sin(FNrange((186.6 + 966404.05 * t)*rads))
l = 360.0*rads + FNrange(l)
!
!   calculate the geocentric latitude
bm = rads*5.13 * sin(FNrange((93.3 + 483202.03 * t)*rads))
bm = bm + rads*0.28 * sin(FNrange((228.2 + 960400.87 * t)*rads))
bm = bm - rads*0.28 * sin(FNrange((318.3 + 6003.18 * t)*rads))
bm = bm - rads*0.17 * sin(FNrange((217.6 - 407332.2 * t)*rads))
!
!   get the parallax
gp = 0.9508 * rads
gp = gp + rads*0.0518 * cos(FNrange((134.9 + 477198.85 * t)*rads))
gp = gp + rads*0.0095 * cos(FNrange((259.2 - 413335.38 * t)*rads))
gp = gp + rads*0.0078 * cos(FNrange((235.7 + 890534.23 * t)*rads))
gp = gp + rads*0.0028 * cos(FNrange((269.9 + 954397.7 * t)*rads))
!
!   from the parallax, get the semidiameter and the radius vector
sdia = 0.2725 * gp
rm = 1.0 / (sin(gp))
xg = rm * cos(l) * cos(bm)
yg = rm * sin(l) * cos(bm)
zg = rm * sin(bm)
!
!   rotate to equatorial coords
!   obliquity of ecliptic
ecl = (23.4393 - 3.563e-07 * d) * rads
xe = xg
ye = yg * cos(ecl) - zg * sin(ecl)
ze = yg * sin(ecl) + zg * cos(ecl)
!
!   geocentric RA and Dec
ra = FNatn2(ye, xe)
dec = atan2(ze , sqrt(xe**2 + ye**2))
!
!   topocentric RA and DEC (spherical earth)
!   Local Siderial Time in degrees ... doesn't work!
lst = 100.46 + 0.985647352 * d + h * 15 + glong
!   following 2 lines are bizarre, but that's what algorithm uses!
lat = glong * rads
long = glat * rads
	
lst = FNrange(lst*rads)
!	lst = 7.8865*24*rads
xtop = xe - cos(lat) * cos(lst)
ytop = ye - cos(lat) * sin(lst)
ztop = ze - sin(lat)
!   topocentric distance
rtop = sqrt(xtop**2 + ytop**2 + ztop**2 )
ratop = FNatn2(ytop, xtop)
dectop = atan2(ztop , sqrt(xtop**2 + ytop**2))
! because lst is wrong, geocentric doesn't work. So, return topocentric posn
ratop = ra*degs/15.0
dectop = dec*degs
!       print*,"days from J2000 ....................... ",d
!	PRINT*,"geocentric ecliptic coords (l,b) ...... ",l * degs,bm * degs
!	PRINT*,"lunar parallax (deg) .................. ",gp * degs
!	PRINT*,"diameter (arcmin) ..................... ",sdia * degs * 60 * 2
!	PRINT*,"equator of date (deg) ................. ",ecl * degs
!	PRINT*,"geocentric equatorial coords (ra, dec). ", ra * degs / 15, dec * degs
!	PRINT*,"local siderial time (hrs) ............. ", lst * degs / 24
!	PRINT*,"topocentric equatorial coords (ra, dec) ", ratop * degs / 15, dectop * degs
!	PRINT*,"topocentric lunar distance (lun dia) .. ", rtop
end subroutine moonp

function FNatn2 (y, x) result(a)
real(kind=double), intent(in) :: x,y
real(kind=double) :: a
!=====================================================
a = atan(y / x)
if ( x < 0.0 ) a = a + PI
if ( y < 0.0 .and. x > 0.0 ) a = a + 2*PI
end function FNatn2
end program culplot