!-------------------------------------------------------------------------
! Copyright 2008, 2009, Roger Hunter
! rogerh@lpbroadband.net
!------------------------------------------------------------------------
! SunGP.exe Version SunGP.exe.1.A.2.2009/01/17
!
! This program is a compiled TrueBasic (TM) adaptation of QBasic
! code written by Keith Burnett
! http://www.stargazing.net/kepler/index.html
!
! and was written in cooperation with E.D.G. for use in his
! earthquake forecasting program,
! ETDPROG Copyright 2004, 2005, 2006, 2007, 2008, 2009, E.D.G.
! seismic@ix.netcom.com
!------------------------------------------------------------------------
! This file is part of ETDPROG.
!
! ETDPROG is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.!
!
! ETDPROG is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with ETDPROG.  If not, see <http://www.gnu.org/licenses/>.
!------------------------------------------------------------------------
PROGRAM sungp (filn$)
!
! ACCEPTS A FILENAME (OR USES DEFAULT FILE) CONTAINING ONE DATE/TIME
DIM res(0:90,2),gr(0:360)

LET pr2$ = "--#.## "
LET pr3$ = "---#.## "
LET pr4$ = "------#.##  "
LET tpi = 2 * pi
LET twopi = tpi
!
!   Get the days to J2000
!   FNday only works between 1901 to 2099 - see Meeus chapter 7
DEF FNday (y, m, d, h, mn)
    LOCAL a,b,c,c1
    LET a = 367 * y
    LET b = int((m + 9)/12)
    LET c = int(7*(y + b) / 4)
    LET c1 = int(275 * m / 9)
    LET e = a - c + c1 + d - 730531.5 + h/24.0 + mn/1440.0
    LET FNday = e
END DEF
!
!   the atn2 function below returns an angle in the range 0 to two pi
!   depending on the signs of x and y.
DEF FNatn2 (y, x)
    IF x = 0 then
       LET a = sgn(y)*rad(90.0)
    ELSE
       LET a = atn(y/x)
    END IF
    IF x < 0 then LET a = a + pi
    IF y < 0 and x > 0 then
       LET a = atn(y/x) + 2*pi
    END IF
    LET FNatn2 = a
END DEF
!
!   the function below returns the true integer part,
!   even for negative numbers
DEF FNipart (x) = SGN(x) * INT(ABS(x))
!
! the function below returns an angle in the range
! 0 to two pi
DEF FNrange (x)
    LET b = x / tpi
    LET a = tpi * (b - FNipart(b))
    IF a < 0 THEN LET a = tpi + a
    LET FNrange = a
END DEF
!
!******* START OF MAIN PROGRAM****************************
!
! OPEN THE OUTPUT FILE

IF filn$ = "" then LET filn$ = "sungpin.dat"
OPEN #2: name filn$, create old

LET ofil$ = "sungpout.dat"
OPEN #3: name ofil$, create newold
SET #3: MARGIN 150
ERASE #3

! INPUT FORMAT "2008/01/01 23:59:10"
INPUT #2: dt$

LET y = val(dt$[1:4])
LET m = val(dt$[6:7])
LET day = val(dt$[9:10])
LET dat$ = dt$[1:10]
!
LET hour = val(dt$[12:13])
LET mins = val(dt$[15:16])
LET sec = val(dt$[18:19])
LET tim$ = dt$[12:19]
!
LET sang = 15*(hour + mins/60 + sec/3600) + 180  ! SUN ANGLE FROM 0 DEGREES
IF sang > 360 then LET sang = sang - 360

LET rlon = abs(glon)
DO while rlon > 90
   LET rlon = rlon - 90
LOOP
LET res(rlon,1) = res(rlon,1) + 1

LET glat = rad(89.9999)
LET glon = rad(0)
!
LET d = FNday (y, m, day, hour, mins)

!-------------------------------------------------
! SUN'S POSITION
!
!   mean longitude of the Sun
LET L = FNrange( rad(280.461 + .9856474*d))
!
!   mean anomaly of the Sun
LET MAs = FNrange( rad(357.528 + .9856003*d))
!
! Find the Earth - Sun distance
LET rs = 1.00014 - .01671 * COS(MAs) - .00014 * COS(2 * MAs)
!
! eccentricity of the sun
LET es = 0.016709 - 1.151e-9*d
!
! argument of perihelion
LET ws = 282.9404 + 4.70935e-5*d
!
!   Ecliptic longitude of the Sun
LET l1 = L + rad(1.915) * SIN(MAs)
LET l2 = rad(.02) * SIN(2 * MAs)
LET lambda = FNrange(L1 + l2)
!
! Obliquity of the ecliptic
LET ecl = rad(23.4393 - 3.563e-7 * d)
!
! Find the RA and DEC of the Sun
LET RAs = FNatn2(COS(ecl)*SIN(lambda), COS(lambda))
LET decs = asin(SIN(ecl)*SIN(lambda))  ! RADIANS
!
! hour angle of Sun
LET LMST = FNrange(rad(280.46061837 + 360.98564736629 * d) + glon)
LET hasun = FNrange(LMST - RAs)
!
!   conversion from hour angle to  Az
LET sinalt = sin(decs)* sin(glat) + cos(decs)*cos(glat)*cos(hasun)
LET y1 = -cos(decs)*cos(glat)*sin(hasun)
LET x1 =  sin(decs) - sin(glat)*sinalt

LET azs = FNatn2(y1, x1)          ! RADIANS
LET azs = azs + pi
LET azs = FNrange(azs)

CALL latlon (decs,azs,glon,lats,lons)  ! DEGREES

!-----------------------------------------------------------------
! Moon positions to a quarter of a degree on the sky
!   From page D76 of Astronomical Almanac
!   
LET t = d / 36525
!
!   calculate the geocentric longitude
LET l = FNrange(rad(218.32 + 481267.883 * t))
LET l = l + rad(6.29 * SIN(FNrange(rad(134.9 + 477198.85 * t))))
LET l = l - rad(1.27 * SIN(FNrange(rad(259.2 - 413335.38 * t))))
LET l = l + rad( .66 * SIN(FNrange(rad(235.7 + 890534.23 * t))))
LET l = l + rad( .21 * SIN(FNrange(rad(269.9 + 954397.7 * t))))
LET l = l - rad( .19 * SIN(FNrange(rad(357.5 + 35999.05 * t))))
LET l = l - rad( .11 * SIN(FNrange(rad(186.6 + 966404.05 * t))))
LET l = FNrange(l)
!
!   calculate the geocentric latitude
LET bm = rad(5.13 * SIN(FNrange(rad(93.3 + 483202.03*t))))
LET bm = bm + rad(.28 * SIN(FNrange(rad(228.2 + 960400.87*t))))
LET bm = bm - rad(.28 * SIN(FNrange(rad(318.3 + 6003.18*t))))
LET bm = bm - rad(.17 * SIN(FNrange(rad(217.6 - 407332.2*t))))
!
! get the parallax
LET gp = rad(.9508)
LET gp = gp + rad(.0518*COS(FNrange(rad(134.9 + 477198.85*t))))
LET gp = gp + rad(.0095*COS(FNrange(rad(259.2 - 413335.38*t))))
LET gp = gp + rad(.0078*COS(FNrange(rad(235.7 + 890534.23*t))))
LET gp = gp + rad(.0028*COS(FNrange(rad(269.9 + 954397.7*t))))
!
! from the parallax, get the semidiameter and the radius vector
LET sdia = .2725 * gp
LET rm = 1 / (SIN(gp))
!
! geocentric ecliptic rectangular coordinates
LET xg = rm * COS(l) * COS(bm)
LET yg = rm * SIN(l) * COS(bm)
LET zg = rm * SIN(bm)
!
! rotate to equatorial coords
LET xe = xg
LET ye = yg * COS(ecl) - zg * SIN(ecl)
LET ze = yg * sin(ecl) + zg * COS(ecl)
!
!   geocentric RA and Dec
LET RAm = FNatn2(ye, xe)          ! RADIANS
LET decm = ATN(ze / SQR(xe*xe + ye*ye))     ! RADIANS
!
! geocentric distance
LET rmg = SQR(xg*xg + yg*yg + zg*zg)   ! EARTH RADII (6378km)
LET rmg = rmg*6378                ! KILOMETERS

!   topocentric RA and DEC (spherical earth)
LET LST = LMST
LET xtop = xe - COS(glat) * COS(lst)
LET ytop = ye - COS(glat) * SIN(lst)
LET ztop = ze - SIN(glat)

LET RAtop = FNatn2(ytop, xtop)    ! RADIANS
LET dectop = ATN(ztop / SQR(xtop*xtop + ytop*ytop))   ! RADIANS
!
LET RAang = deg(RAs-RAtop)
IF RAang <= -180 then LET RAang = RAang + 360
IF RAang > 180 then LET RAang = RAang - 360

! topocentric distance
LET rmtop = SQR(xtop*xtop + ytop*ytop + ztop*ztop)    ! EARTH RADII (6378km)
LET rmtop = rmtop*6378            ! KILOMETERS

LET LHA = LST - RAm               ! radians

! Compute essential sines, cosines and tangents
LET Sin_LHA = Sin((LHA))
LET Cos_LHA = Cos((LHA))

LET Sin_B = Sin(gLat)
LET Cos_B = Cos(gLat)

LET Sin_Decl = Sin(Dectop)
LET Cos_Decl = Cos(Dectop)
LET Tan_Decl = Tan(Dectop)

! Compute horizon azimuth reckoned clockwise from south.
LET Azm = FNatn2(Sin_LHA, Cos_LHA*Sin_B - Tan_Decl*Cos_B)
LET azm = FNrange(azm)            ! convert to north (radians)

CALL latlon (dectop,azm,glon,latm,lonm)     ! DEGREES

LET azang = deg(azs - azm)
IF azang <= -180 then LET azang = azang + 360
IF azang > 180 then LET azang = azang - 360

! PRINT RESULTS IN DECIMAL FORM
LET mag$ = str$(mag)
IF len(mag$) = 1 then LET mag$ = mag$ & ".0"

LET lats = round(lats,2)
LET lons = round(lons,2)

LET latm = round(latm,2)
LET lonm = round(lonm,2)

! GP LAT   
CALL gp2 (lats,lons,latm,lonm,Y)

LET redl = lonm
DO while redl > 90
   LET redl = redl - 90
LOOP

! GP LON, REAL
CALL gp1 (lats,lons,latm,lonm,a1,a2)
LET ndx = int(a2)                   ! 0-360
LET gr(ndx) = gr(ndx) + 1
LET ra = abs(a1)                    ! 0-90
IF ra > 90 then LET ra = ra - 90

LET res(ra,2) = res(ra,2) + 1

LET rb = a2
DO while rb > 90
   LET rb = rb - 90
LOOP

LET azs = deg(azs)
IF azs < 0 then LET azs = azs + 360
IF azs >= 360 then LET azs = azs - 360
LET azs = round(azs,2)

LET azm = deg(azm)
IF azm < 0 then LET azm = azm + 360
IF azm >= 360 then LET azm = azm - 360
LET azm = round(azm,2)
!
LET azang = round(azang,2)

! DISTANCE TO SUN
LET rs = round(rs,3)

! DISTANCE TO MOON
LET rmg = round(rmg,0)

! GP STRENGTH
CALL gp3 (lats,lons,latm,lonm,H)

PRINT #3: dat$ & "  UTC format input date"
PRINT #3: "  " & tim$ & "  UTC format input time" 
PRINT #3, USING pr4$: lats;
print #3: "Subsolar point latitude"

PRINT #3, USING pr4$: lons;
print #3: "Subsolar point longitude"

PRINT #3, USING pr4$: latm;
print #3: "Sublunar Point latitude"

PRINT #3, USING pr4$: lonm;
print #3: "Sublunar Point longitude"

PRINT #3, using pr4$: Y;
print #3: "True Gravity Point latitude"

print #3, using pr4$: a2;
print #3: "True Gravity Point longitude"

PRINT #3, using pr4$: redl;
print #3: "0 to 89 degree range value for Sublunar Point longitude"

PRINT #3, using pr4$: rb;
print #3: "0 to 89 degree range value for Gravity Point longitude"

PRINT #3, using pr4$: azang;
print #3: "Sun - Earth - Moon angle"

print #3: "  ";
PRINT #3, using "---#.###  ": rs;
print #3: "Distance between Earth and sun (astronomical units)"

PRINT #3, using "##########  ": rmg;
print #3: "Distance between Earth and moon (kilometers)"

PRINT #3, using pr4$: H;
print #3: "Combined sun and moon gravities at the Gravity Point"

CLOSE #2
CLOSE #3

!**************************************
SUB latlon (dec,az,glon,lat,lon)
    LET lat = deg(dec)
    LET lon = FNrange(az - glon)
    LET lon = deg(lon)

    IF lon < 0 then LET lon = lon + 360
    IF lon > 360 then
       LET lon = lon - 360
    END IF

END SUB

!**************************************
SUB gp1 (al,am,an,ao,a2,a3)
    ! GRAVITY POINT LONGITUDES

    LOCAL ap,aq,ar,as,at,au,aw,ax,bj,bk,bn,bo

    LET ap = rad(al)
    LET aq = rad(am)

    LET ar = rad(an)
    LET as = rad(ao)

    LET at = cos(ap)*sin(aq)
    LET au = cos(ap)*cos(aq)

    LET aw = cos(ar)*cos(as)
    LET ax = cos(ar)*sin (as)

    LET bj = au + aw*2.5
    IF bj = 0 then LET bj = bj + 0.000001

    LET bk = at + ax*2.5

    LET bn = atn(bk/bj)
    LET bn = deg(bn)

    LET a3 = bn                   ! REAL 0 - 360
    IF bj < 0 then LET a3 = a3 + 180
    IF (bj > 0 and bk < 0) then LET a3 = a3 + 360
    IF a3 >= 360 then LET a3 = a3 - 360
    IF a3 < 0 then LET a3 = a3 + 360

    ! NEXT FEW LINES CHANGED
    LET a2 = a3                   ! USE  0 - 90
    DO while a2 >= 90
       LET a2 = a2 - 90
    LOOP

END SUB

!**********************************************
SUB gp2 (H,I,J,K,Y)
    LOCAL L,M,N,O,P,Q,R,S,T,U,V,W,X

    LET L = rad(H)
    LET M = rad(I)

    LET N = rad(J)
    LET O = rad(K)

    LET P = cos(L)*cos(M)
    LET Q = cos(L)*sin(M)
    LET R = sin(L)

    LET S = cos(N)*cos(O)
    LET T = cos(N)*sin(O)

    LET U = sin(N)

    LET V = P + S*2.5
    LET W = Q + T*2.5
    LET X = R + U*2.5

    LET Y = atn(X/sqr(V^2 + W^2))
    LET Y = deg(Y)

END SUB

!***********************************************
SUB gp3 (C,D,E,F,H)
    ! GRAVITY POINT STRENGTH
    LOCAL I,Ia,Ib,J,Ja,Jb,K,Ka,N,O

    IF (d-f) > 180 then
       LET Ia = 360
    ELSE
       LET Ia = 0
    END IF
    IF (f-d) > 180 then
       LET Ib = 360
    ELSE
       LET Ib = 0
    END IF
    LET I = ABS(D-F - Ia) + Ib

    LET Ja = 0
    IF (I <= 90) then
       LET Ja = 1/5400
    END IF
    LET J = 0.38*ABS(C-E)*(90-I)*Ja

    LET Ka = 0
    IF (I > 90) then
       LET Ka = 1/5400
    END IF

    LET K = 1.32*ABS((C+E)*(I-90))*Ka
    LET N = 2.5*cos(rad(F)) + cos(rad(D))
    LET O = 2.5*sin(rad(F)) + sin(rad(D))

    LET H = SQR(N^2 + O^2) + J + K

END SUB



END