Deltares / hatyan

Harmonic tidal analysis and prediction
https://deltares.github.io/hatyan/
GNU General Public License v3.0
13 stars 2 forks source link

`astrog_moonriseset` fails for some longitudes #323

Open veenstrajelmer opened 3 months ago

veenstrajelmer commented 3 months ago
import pytz
import datetime as dt
import pandas as pd
import matplotlib.pyplot as plt
plt.close('all')
import hatyan

# script settings
timeStart = dt.datetime(2022,1,20)
timeEnd   = dt.datetime(2022,1,27)
dT_fortran = False #True is best comparison to fortran, False is more precise 
tz_GMT = 'UTC' # UTC/GMT timezone
try: #only works with pandas 1.2.0 or higher, but 1.1.5 is highest available pandas version in Python 3.6.12 (which is the highest available Python in RHEL)
    tz_MET = 'UTC+01:00' #for timezone conversion to UTC+01:00 (Dutch wintertime)
    pytz.timezone(tz_MET)
except:
    tz_MET = dt.timezone(dt.timedelta(hours=1)) #for timezone conversion to UTC+01:00 (Dutch wintertime)
tz_EurAms = 'Europe/Amsterdam' # for conversion to local timezone, including daylight saving time (DST)

pdtocsv_kwargs = dict(index=False, sep=';', date_format='%Y%m%d %H%M%S', float_format='%9.5f', na_rep='--:--')

# sunrise and -set
# ook uitgegeven door KNMI voor lon=5, lat=52 (ronde coordinaten, rotonde bij IJsselstein) voor tabel met zon op/onder en begin seizoenen:
# https://cdn.knmi.nl/ckeditor_assets/attachments/170/Tijden_van_zonopkomst_en_-ondergang_2022.pdf

print('sun')
for lon in []:#np.arange(-180,180+1,45): #30 degrees is 30/360*24=2 hours
    print()
    print(lon)
    tz_LOCAL = dt.timezone(dt.timedelta(hours=lon/360*24))
    sunriseset_python = hatyan.astrog_sunriseset(tFirst=timeStart, tLast=timeEnd, dT_fortran=dT_fortran, tzone=tz_GMT, lon=lon, lat=52.0)
    print(sunriseset_python.iloc[-2:])

lat,lon =  54.7165, 135.3084 #Vladivostok #this one crashes for longer time periods (rates of increase go off track). Sort of solved by switching signs of ALTMOO RATE in astrac, but moonrise/set are then switched and it crashes for lat,lon=50,45
lat,lon =  55, 45.3084 #fake
lat,lon = -33.8688, 151.2093 #sydney
#lat,lon =  52.1561,   5.3878 #amersfoort
#lat,lon =  51.47869,  -0.01080 #greenwich
sunriseset_python = hatyan.astrog_moonriseset(tFirst=timeStart, tLast=timeEnd, dT_fortran=dT_fortran, tzone='UTC+10:00', lon=lon,lat=lat)
#sunriseset_python = hatyan.astrog_sunriseset(tFirst=timeStart, tLast=timeEnd, dT_fortran=dT_fortran, tzone='UTC', lon=lon,lat=lat)
datesmoonrise = pd.DatetimeIndex(sunriseset_python.loc[:,'datetime'].dt.tz_localize(None))
print(sunriseset_python)
astrabOutput = hatyan.astrog.astrab(datesmoonrise,dT_fortran=dT_fortran,lon=lon,lat=lat)
#print(astrabOutput['ALTSUN'])
#print(astrabOutput['ALTMOO'])
#print(astrabOutput['ALTMOO'])

print('moon')
for lon in range(-180,180+1,45): #30 degrees is 30/360*24=2 hours suntime
    print()
    print(lon)
    tz_LOCAL = dt.timezone(dt.timedelta(hours=lon/360*24))
    try:
        lat=50
        tzone_local = dt.timezone(dt.timedelta(hours=lon/360*24))
        sunriseset_python = hatyan.astrog_moonriseset(tFirst=timeStart, tLast=timeEnd, dT_fortran=dT_fortran, tzone=tz_GMT, lon=lon, lat=lat) 
        #print('SUCCES')
        print(sunriseset_python.iloc[:3])
        datesmoonrise = pd.DatetimeIndex(sunriseset_python.loc[:,'datetime'].dt.tz_localize(None))
        astrabOutput = hatyan.astrog.astrab(datesmoonrise,dT_fortran=dT_fortran,lon=lon,lat=lat)
        print(astrabOutput['ALTMOO'])
    except Exception as e:
        print(f'FAILED: {e}')

Gives:

sun
                    datetime  type  type_str
0  2022-01-20 06:43:35+10:00     1  moonrise
1  2022-01-20 20:52:34+10:00     2   moonset
2  2022-01-21 07:44:24+10:00     1  moonrise
3  2022-01-21 21:23:46+10:00     2   moonset
4  2022-01-22 08:45:10+10:00     1  moonrise
5  2022-01-22 21:52:49+10:00     2   moonset
6  2022-01-23 09:46:03+10:00     1  moonrise
7  2022-01-23 22:21:01+10:00     2   moonset
8  2022-01-24 10:47:46+10:00     1  moonrise
9  2022-01-24 22:49:46+10:00     2   moonset
10 2022-01-25 11:51:16+10:00     1  moonrise
11 2022-01-25 23:20:36+10:00     2   moonset
12 2022-01-26 12:57:30+10:00     1  moonrise
13 2022-01-26 23:55:18+10:00     2   moonset
moon

-180
                   datetime  type  type_str
0 2022-01-20 06:34:02+00:00     2   moonset
1 2022-01-20 21:42:34+00:00     1  moonrise
2 2022-01-21 07:47:33+00:00     2   moonset
[-0.9754931  -0.65127786 -0.98680884 -0.64586638 -0.99779241 -0.64229981
 -1.00121652 -0.64429865 -1.00294253 -0.65300772 -1.00023903 -0.66385949
 -0.99109837 -0.67917   ]

-135
FAILED: Iteration step resulted in time changes larger than 24 hours (max 59.41 hours), try using a lower latitude

-90
                   datetime  type  type_str
0 2022-01-20 00:15:47+00:00     1  moonrise
1 2022-01-20 15:37:39+00:00     2   moonset
2 2022-01-21 01:29:08+00:00     1  moonrise
[-0.97111945 -0.65345865 -0.98361091 -0.64681695 -0.99440636 -0.64425597
 -1.00054758 -0.64431841 -1.00407278 -0.6488786  -1.00008469 -0.65981283
 -0.99351902 -0.67436023]

-45
                   datetime  type  type_str
0 2022-01-20 12:35:07+00:00     2   moonset
1 2022-01-20 22:19:56+00:00     1  moonrise
2 2022-01-21 12:54:02+00:00     2   moonset
[-0.65597334 -0.98178981 -0.64837047 -0.99285961 -0.64416026 -1.00093252
 -0.64387019 -1.0045024  -0.64758008 -1.00179117 -0.65793708 -0.99421506
 -0.67437941]

0
                   datetime  type  type_str
0 2022-01-20 09:32:30+00:00     1  moonrise
1 2022-01-20 19:10:44+00:00     2   moonset
2 2022-01-21 09:51:49+00:00     1  moonrise
[-0.65533445 -0.98078679 -0.64784179 -0.99203119 -0.64386654 -1.00136471
 -0.64276793 -1.0026784  -0.64745481 -1.001186   -0.65789562 -0.99500145
 -0.67022731]

45
FAILED: Iteration step resulted in time changes larger than 24 hours (max 76.85 hours), try using a lower latitude

90
                   datetime  type  type_str
0 2022-01-20 03:27:06+00:00     2   moonset
1 2022-01-20 12:52:22+00:00     1  moonrise
2 2022-01-21 03:47:17+00:00     2   moonset
[-0.65847258 -0.9773919  -0.64981116 -0.98965461 -0.64423219 -0.99865603
 -0.64282382 -1.00353081 -0.64736334 -1.00200136 -0.65575052 -0.99770357
 -0.66691263 -0.98879438]

135
FAILED: Iteration step resulted in time changes larger than 24 hours (max 28.29 hours), try using a lower latitude

180
                   datetime  type  type_str
0 2022-01-20 06:34:02+00:00     2   moonset
1 2022-01-20 21:42:34+00:00     1  moonrise
2 2022-01-21 07:47:33+00:00     2   moonset
[-0.9754931  -0.65127786 -0.98680884 -0.64586638 -0.99779241 -0.64229981
 -1.00121652 -0.64429865 -1.00294253 -0.65300772 -1.00023903 -0.66385949
 -0.99109837 -0.67917   ]

Probably also for hatyan.astrog_sunriseset()