# NYCSpiral.py
# Andrew Davison, ad@coe.psu.ac.th, July 2025
# input: 80
'''
Draw a long Rhumb line starting from NYC, with the bearing
supplied by the user. It will spiral towards a pole.
'''

import cartopy.crs as ccrs
import cartopy.feature as cf
import matplotlib.pyplot as plt
import nav

MAX_COORDS = 20000

def plotCity(ax, city, lat, lon):
  # Add marker
  ax.plot(lon, lat, 'o', color='green', 
         transform=ccrs.PlateCarree())
  # label city
  plt.text(lon, lat, " "+city, color='green', 
         transform=ccrs.PlateCarree())


def rhumbbearingCoords(lat1, lon1, bearing, dist=100): # 100
  nCoords = 0
  coords = [(lat1, lon1)]
  currLat, currLon = lat1, lon1
  while (abs(currLat) < 89) and (nCoords < MAX_COORDS):
    destLat, destLon =  \
        nav.rhumbDestCoord(currLat, currLon, dist, bearing)
    coords.append((destLat, destLon))
    currLat, currLon = destLat, destLon
    nCoords += 1
  return coords


# get coordinates for NYC
bearing = float(input("Bearing from NYC? "))
startCity = 'New York'
cities = nav.loadCityData()
lat1, lon1 = nav.readCoord(cities, startCity)
print(f"{startCity} (lat,lon): ({lat1:.2f}, {lon1:.2f}); {bearing:.2f} degs")

plt.figure("NYC Spiral", figsize=(6,6))

proj = ccrs.Orthographic(lon1, 40)
proj.threshold = proj.threshold/20  # finer threshold
ax = plt.axes(projection = proj)


ax.gridlines(draw_labels=True, x_inline=True, ls='--',
             xlocs=range(-180, 181, 30),
             ylocs=range(-90, 91, 10))
ax.coastlines()
ax.set_global()  # show the entire globe
ax.add_feature(cf.BORDERS, ls=':')


# draw the Rhumb line trip as a series of coordinates
coords = rhumbbearingCoords(lat1, lon1, bearing)
'''
print("Rhumb Line (lat,lon)s:")
for cLat, cLon in coords:
  print(f"  ({cLat:.2f}, {cLon:.2f})", end='')
print()
'''
lats, lons = zip(*coords)
plt.plot(lons, lats, color='red',
                      transform=ccrs.Geodetic())
       # to deal with curvature between points

plotCity(ax, startCity, lat1, lon1)

plt.title("Rhumb line, bearing= " + str(bearing) + "\n")
plt.show()