CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

| Download

Calculations and F2 height data pulls to determine the radiation angle of the Project TouCans dipole.

Project: Catenary
Views: 45
Image: ubuntu2204
Tweet
Share
Share
import numpy
import math
import requests
import datetime

deg2rad = math.pi/180
rad2deg = 180/math.pi

def cartesian_x(f,l):
    #f = latitude, l = longitude
    return (math.cos(f*deg2rad)*math.cos(l*deg2rad))

def cartesian_y(f,l):
    #f = latitude, l = longitude
    return (math.cos(f*deg2rad)*math.sin(l*deg2rad))

def cartesian_z(f,l):
    #f = latitude, l = longitude
    return (math.sin(f*deg2rad))

def cross_x(x, y, z, i,j,k):
    return ((y*k)-(z*j))
def cross_y(x, y, z, i,j,k):
    return ((z*i)-(x*k))
def cross_z(x, y, z, i,j,k):
    return ((x*j)-(y*i))

def spherical_lat(x,y,z):
    r = math.sqrt(x*x + y*y)
    #Omitting the special cases because points will always
    #be separated for this application
    return (math.atan2(z, r)*rad2deg) # return degrees

def spherical_lng(x,y,z):
    #Omitting the special cases because points will always
    #be separated for this application
    return (math.atan2(y, x)*rad2deg) # return degrees


#Let's make a sphere for the Earth

earth = sphere(color='pink', center=(0,0,0), size=1, opacity=0.4)
earth

#Now let's get the height of the nearby ionosphere so that we can make another sphere
#First, the radius of the Earth in km. We'll need this
earth_rad_km = 6371

#What was the height of the f2 when the spot happened?
fromDate = "2023-08-21T03:43:00"
endDate =  "2023-08-21T05:43:00"
start_date_win = fromDate#.strftime("%Y.%m.%d %H:%M:%S")
start_date_win = start_date_win.replace(" ", "%20")
    #print(start_date_win)
end_date_win = endDate#.strftime("%Y.%m.%d %H:%M:%S")
end_date_win = end_date_win.replace(" ", "%20")
    #print("start_date_win " + start_date_win)
print("end_date_win " + end_date_win)
#https://github.com/hcarter333/rm-rbn-history/blob/0e4ed3b1867917be88a9122839c91b781d2fb54b/ionodata.py#L16
iono_url = "https://lgdc.uml.edu/common/DIDBGetValues?ursiCode=PA836&charName=MUFD,hF2,hmF2&fromDate=" + start_date_win + "&toDate=" + end_date_win

#iono_url = "https://lgdc.uml.edu/common/DIDBGetValues?ursiCode=PA836&charName=MUFD,hmF2&fromDate=2024.01.27T11:00&toDate=2024.01.27T14:23"

iodata=requests.get(iono_url)
#print(iodata.text)
iono_data = iodata.text
iono_lines = iono_data.split("\n")
data_found = 0
for data_line in iono_lines:
    print(data_line)
#    if(data_found == 1):
#        #We've arrived at the data
#        #Split on space and print the max F2
#        data_fields = data_line.split()
#        hmF2 = data_fields[4]
#        print("hmF2 " + data_fields[4])
        #break
#we're going to use the closest time to the beginning of the range for now
#    if(data_line.find("#Time") != -1):
#        data_found = 1

#show(earth)



#This might not be the most efficient way, but let's try to
#get the 18 degree more accurately
#lat, lng for tx
t_lat = 37.7248952200944
t_lng = -122.422936174405
#lat,lng for landing point
r_lat = 55.24556451189041
r_lng = -110.29305571520665

tx_x = cartesian_x(t_lat,t_lng)
tx_y = cartesian_y(t_lat,t_lng)
tx_z = cartesian_z(t_lat,t_lng)
rx_x = cartesian_x(r_lat,r_lng)
rx_y = cartesian_y(r_lat,r_lng)
rx_z = cartesian_z(r_lat,r_lng)

#Take the cross product of the two positions
g_x = cross_x(tx_x, tx_y, tx_z, rx_x, rx_y, rx_z)
g_y = cross_y(tx_x, tx_y, tx_z, rx_x, rx_y, rx_z)
g_z = cross_z(tx_x, tx_y, tx_z, rx_x, rx_y, rx_z)

#Now, get the magnitude of the cross product
mag_g = math.sqrt(g_x**2+g_y**2+g_z**2)
mag_tx = math.sqrt(tx_x**2+tx_y**2+tx_z**2)
print("magnitude of cross product" + str(mag_g))
print("magnitude of cross product" + str(mag_tx))
theta = math.asin(mag_g)*rad2deg
print("swept_angle = " + str(theta))

tx_cos = math.cos(deg2rad*theta)
F2_min = earth_rad_km*((1-tx_cos)/tx_cos)
print("Calculated minimum height of F2", str(F2_min))
print("radius of earth in km ", earth_rad_km)

print("Now let's look at the F2 height at Idaho")
iono_url = "https://lgdc.uml.edu/common/DIDBGetValues?ursiCode=IF843&charName=MUFD,hF2,hmF2&fromDate=" + start_date_win + "&toDate=" + end_date_win

#iono_url = "https://lgdc.uml.edu/common/DIDBGetValues?ursiCode=PA836&charName=MUFD,hmF2&fromDate=2024.01.27T11:00&toDate=2024.01.27T14:23"

iodata=requests.get(iono_url)
#print(iodata.text)
iono_data = iodata.text
iono_lines = iono_data.split("\n")
data_found = 0
for data_line in iono_lines:
    print(data_line)

print("\n\n\nAnd then from Thule")
iono_url = "https://lgdc.uml.edu/common/DIDBGetValues?ursiCode=THJ76&charName=MUFD,hF2,hmF2&fromDate=" + start_date_win + "&toDate=" + end_date_win

#iono_url = "https://lgdc.uml.edu/common/DIDBGetValues?ursiCode=PA836&charName=MUFD,hmF2&fromDate=2024.01.27T11:00&toDate=2024.01.27T14:23"

iodata=requests.get(iono_url)
#print(iodata.text)
iono_data = iodata.text
iono_lines = iono_data.split("\n")
data_found = 0
for data_line in iono_lines:
    print(data_line)
/ext/sage/10.3/local/var/lib/sage/venv-python3.11.1/lib/python3.11/site-packages/scikits/__init__.py:1: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html __import__("pkg_resources").declare_namespace(__name__) /ext/sage/10.3/local/var/lib/sage/venv-python3.11.1/lib/python3.11/site-packages/scikits/__init__.py:1: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('scikits')`. Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages __import__("pkg_resources").declare_namespace(__name__)
3D rendering not yet implemented
end_date_win 2023-08-21T05:43:00 # Global Ionospheric Radio Observatory # GIRO Tabulated Ionospheric Characteristics, Version 1.0 Revision B # Generated by DIDBGetValues on 2024-03-27T18:16:01.514Z # # Location: GEO 34.8N 239.5E, URSI-Code PA836 PT ARGUELLO # Instrument: Ionosonde, Model: DPS-4D # # Query for measurement intervals of time: # 2023-08-21T03:43:00.000Z - 2023-08-21T05:43:00.000Z # # Data Selection: # CS is Autoscaling Confidence Score (from 0 to 100, 999 if manual scaling, -1 if unknown) # MUFD [MHz] - Maximum usable frequency for ground distance D # hF2 [km] - Minimum virtual height of F2 trace # hmF2 [km] - Peak height F2-layer # Distance D for MUF calculations: 3000 km # # All GIRO measurements are released under CC-BY-NC-SA 4.0 license # Please follow the Lowell GIRO Data Center RULES OF THE ROAD # https://ulcar.uml.edu/DIDBase/RulesOfTheRoadForDIDBase.htm # Requires acknowledgement of PA836 data provider # #Time CS MUFD QD hF2 QD hmF2 QD 2023-08-21T03:45:00.000Z 0 19.69 // 247.5 // --- __ 2023-08-21T03:52:30.000Z 80 19.38 // 247.5 // 310.9 // 2023-08-21T04:00:00.000Z 80 19.40 // 240.0 // 312.6 // 2023-08-21T04:07:30.000Z 80 18.94 // 240.0 // 309.5 // 2023-08-21T04:15:00.000Z 0 18.94 // 245.0 // --- __ 2023-08-21T04:22:30.000Z 35 18.96 // 247.5 // 282.5 // 2023-08-21T04:30:00.000Z 0 18.24 // 242.5 // --- __ 2023-08-21T04:37:30.000Z 0 17.49 // 265.0 // --- __ 2023-08-21T04:45:00.000Z 0 17.43 // 255.0 // --- __ 2023-08-21T04:52:30.000Z 0 16.83 // 232.5 // --- __ 2023-08-21T05:00:00.000Z 0 16.85 // 240.0 // --- __ 2023-08-21T05:07:30.000Z 80 16.68 // 245.0 // 335.5 // 2023-08-21T05:15:00.000Z 60 15.96 // 252.5 // 277.9 // 2023-08-21T05:22:30.000Z 80 16.10 // 250.0 // 330.3 // 2023-08-21T05:30:00.000Z 0 15.71 // 255.0 // --- __ 2023-08-21T05:37:30.000Z 60 15.00 // 257.5 // 307.9 // magnitude of cross product0.3312538812032753 magnitude of cross product1.0 swept_angle = 19.34489864958923 Calculated minimum height of F2 381.2190729008156 radius of earth in km 6371 Now let's look at the F2 height at Idaho # Global Ionospheric Radio Observatory # GIRO Tabulated Ionospheric Characteristics, Version 1.0 Revision B # Generated by DIDBGetValues on 2024-03-27T18:16:01.817Z # # Location: GEO 43.81N 247.32E, URSI-Code IF843 IDAHO NATIONAL LAB # Instrument: Ionosonde, Model: DPS-4D # # Query for measurement intervals of time: # 2023-08-21T03:43:00.000Z - 2023-08-21T05:43:00.000Z # # Data Selection: # CS is Autoscaling Confidence Score (from 0 to 100, 999 if manual scaling, -1 if unknown) # MUFD [MHz] - Maximum usable frequency for ground distance D # hF2 [km] - Minimum virtual height of F2 trace # hmF2 [km] - Peak height F2-layer # Distance D for MUF calculations: 3000 km # # All GIRO measurements are released under CC-BY-NC-SA 4.0 license # Please follow the Lowell GIRO Data Center RULES OF THE ROAD # https://ulcar.uml.edu/DIDBase/RulesOfTheRoadForDIDBase.htm # Requires acknowledgement of IF843 data provider # #Time CS MUFD QD hF2 QD hmF2 QD 2023-08-21T03:45:00.000Z 70 20.11 // 237.0 // 296.1 // 2023-08-21T04:00:00.000Z 45 19.95 // 240.0 // 275.6 // 2023-08-21T04:15:00.000Z 70 20.68 // 236.5 // 303.4 // 2023-08-21T04:30:00.000Z 70 19.77 // 240.0 // 291.6 // 2023-08-21T04:45:00.000Z 40 17.59 // 237.0 // 291.4 // 2023-08-21T05:00:00.000Z 70 16.59 // 265.0 // 335.7 // 2023-08-21T05:15:00.000Z 70 16.47 // 253.0 // 325.3 // 2023-08-21T05:30:00.000Z 70 15.99 // 247.0 // 313.3 // And then from Thule # Global Ionospheric Radio Observatory # GIRO Tabulated Ionospheric Characteristics, Version 1.0 Revision B # Generated by DIDBGetValues on 2024-03-27T18:16:02.063Z # # Location: GEO 76.54N 291.56E, URSI-Code THJ76 THULE # Instrument: Ionosonde, Model: DPS-4D # # Query for measurement intervals of time: # 2023-08-21T03:43:00.000Z - 2023-08-21T05:43:00.000Z # # Data Selection: # CS is Autoscaling Confidence Score (from 0 to 100, 999 if manual scaling, -1 if unknown) # MUFD [MHz] - Maximum usable frequency for ground distance D # hF2 [km] - Minimum virtual height of F2 trace # hmF2 [km] - Peak height F2-layer # Distance D for MUF calculations: 3000 km # # All GIRO measurements are released under CC-BY-NC-SA 4.0 license # Please follow the Lowell GIRO Data Center RULES OF THE ROAD # https://ulcar.uml.edu/DIDBase/RulesOfTheRoadForDIDBase.htm # Requires acknowledgement of THJ76 data provider # #Time CS MUFD QD hF2 QD hmF2 QD 2023-08-21T03:45:00.000Z 60 15.25 // 210.0 // 226.6 // 2023-08-21T03:52:30.000Z 40 14.09 // 217.5 // 271.1 // 2023-08-21T04:00:00.000Z 60 14.01 // 230.0 // 271.3 // 2023-08-21T04:07:30.000Z 20 13.32 // 222.5 // 215.1 // 2023-08-21T04:15:00.000Z 60 13.35 // 227.5 // 324.4 // 2023-08-21T04:22:30.000Z 70 13.05 // 225.0 // 302.3 // 2023-08-21T04:30:00.000Z 20 13.13 // 318.6 // 251.9 // 2023-08-21T04:37:30.000Z 60 14.83 // 242.5 // 323.3 // 2023-08-21T04:45:00.000Z 60 16.93 // 280.0 // 289.2 // 2023-08-21T04:52:30.000Z 60 17.15 // 252.5 // 299.9 // 2023-08-21T05:00:00.000Z 60 16.92 // 232.5 // 284.3 // 2023-08-21T05:07:30.000Z 60 17.16 // 222.5 // 257.2 // 2023-08-21T05:15:00.000Z 70 16.35 // 257.5 // 267.8 // 2023-08-21T05:22:30.000Z 70 15.26 // 272.5 // 293.7 // 2023-08-21T05:30:00.000Z 70 13.45 // 227.5 // 319.0 // 2023-08-21T05:37:30.000Z 70 13.99 // 292.5 // 302.6 // 




#catenary(a,x) = a*cosh(x/a)
#a2 = plot(catenary(0.2,x), 0, 0.5, color='red', legend_label='a=0.2', legend_color='red')
#a3 = plot(catenary(0.3,x), 0, 0.5, color='green', legend_label='a=0.3', legend_color='green')
#a4 = plot(catenary(0.4,x), 0, 0.5, color='blue', le
exam(x)=-2*(x^4)-(5*(x^3))-(x^2)-(8*x)-6
hamie=plot(exam(x), -3, 3, color='green', ymin=-30,ymax=40)
hamie