1039.sagews

#1st order differential equations (o.d.e) for refraction:
#dr/ds = u, where u = cos(theta)
#du/ds = (1-u^2)*(1/r+1/n*dn/dr)
#n(r) = 1+10^-6*a*exp(-b*r), where n is index of refraction (Exp Model)
#dn/dr = 1+10^-6*a*-b*exp(-b*r)

#Pick initial condition:
#theta = 30 degrees, in radians: pi/6
#u = cos(pi/6) = 0.86602
#r = 1.7(per m)*10^4

radian = var('radian')
radian = pi/6
u = cos(radian).n(digits = 5)
u #Returns initial condition of u

0.86602

#From problem 3, where r = h:
#a = 340
#b = 1.4090(per m)*10^4

#Use 2nd Order Runge-Kutta to integrate equation du/ds = (1-u^2)*(1/r+1/n*dn/dr)

def udot(r,u):
    a = 340
    b = 1.4090
    return (1-(u)^2)*((1)/(r))+((1)/(1+10^-6*a*exp(-b*r)))*(1+10^-6*a*-b*exp(-b*r))

import scipy #Transfer files from scipy (open source library of algorithms and mathematical tools from Python programming language) into this worksheet.

def rk2(u0,r0,rn,rt): #2nd Order Runge-Kutta defined.

    radial_dis = scipy.arange(r0,rn,rt) #Creates array range (arange) of radial_dis values.
    u=[] #Stores u = dr/ds or derivative of radial distance with respect to (s) measurement along               
         #the ray path.
    for r in radial_dis:
        u.append(u0) #Attach results (u0) to variable u.
        k1=udot(r,u0)
        k2=udot(r+0.5*rt,u0+0.5*k1)
        u0=u0+(k1+k2)/2
    return radial_dis,u

rk2(0.86602, 1.71, 3, 0.1)

(array([ 1.71, 1.81, 1.91, 2.01, 2.11, 2.21, 2.31, 2.41, 2.51, 2.61, 2.71, 2.81, 2.91]), [0.866020000000000, 1.6347991621959448, 1.6893601550920172, 1.7183588905432456, 1.7466747042757742, 1.7745479694167945, 1.8020015563352152, 1.829053772212458, 1.855721580897002, 1.8820207910761586, 1.9079661643731738, 1.9335715093204211, 1.9588497647662324])

r,u = rk2(0.86602, 1.71, 3, 0.1)

import numpy as np #Transfer files from fundamental library numpy (N-dimensional array manipulation) used with Python.
import matplotlib.pyplot as plt #Transfer files from matplotlib.pyplot (plotting library) used with Python language. 
#Index of Refraction is always greater than 1.  Index of Refraction approaches 1 from above.

x = u
y = r

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, y)
ax.set_xlim(0.76602, 2.2)
ax.set_ylim(1.6, 3.5)
ax.set_xlabel('Angle of Incidence (u) in radians')
ax.set_ylabel('Radial Distance (r) in (m)*10^4 above mean sea level')
ax.set_title('Numerically Simulated Differential Equations of Refraction')
ax.grid()

plt.savefig('myfig.png')

#From the graph above as the measurement along the ray path (s) changes, the radial distance (r) and the angle of incidence (u) changes.

#Use 4th Order Runge-Kutta to integrate equation du/ds = (1-u^2)*(1/r+1/n*dn/dr)

import scipy #Transfer files from scipy (open source library of algorithms and mathematical tools from Python programming language) into this worksheet.

def rk4(u0,r0,rn,rt): #4th Order Runge-Kutta defined.

    radial_dis = scipy.arange(r0,rn,rt) #Creates array range (arange) of radial_dis values.
    u=[] #Stores u = dr/ds or derivative of radial distance with respect to (s) measurement along               
         #the ray path.

    for r in radial_dis:
   	u.append(u0) #Attach results (x0) at end of variable x.
   	k1 = udot(r,u0)
   	k2 = udot(r+0.5*rt, u0 + 0.5*rt*k1)
   	k3 = udot(r+0.5*rt, u0 + 0.5*rt*k2)
   	k4 = udot(r+rt, u0 + rt*k3)
   	u0 = u0+(rt/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4)
    return radial_dis,u

rk4(0.86602, 1.71, 3, 0.1)

(array([ 1.71, 1.81, 1.91, 2.01, 2.11, 2.21, 2.31, 2.41, 2.51, 2.61, 2.71, 2.81, 2.91]), [0.866020000000000, 0.97460443773712235, 1.0719728174147416, 1.1593986649779751, 1.2381088562350016, 1.3092357190759094, 1.3737956206639728, 1.4326839736496095, 1.486679711579751, 1.5364546840656059, 1.5825851336925272, 1.625563573562677, 1.6658101345483927])

r1,u1 = rk4(0.86602, 1.71, 3, 0.1)

import numpy as np #Transfer files from fundamental library numpy (N-dimensional array manipulation) used with Python.
import matplotlib.pyplot as plt #Transfer files from matplotlib.pyplot (plotting library) used with Python language. 
#Index of Refraction is always greater than 1.  Index of Refraction approaches 1 from above.

x = u1
y = r1

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, y)
ax.set_xlim(0.76602, 1.8)
ax.set_ylim(1.6, 3.5)
ax.set_xlabel('Angle of Incidence (u) in radians')
ax.set_ylabel('Radial Distance (r) in (m)*10^4 above mean sea level')
ax.set_title('Numerically Simulated Differential Equations of Refraction')
ax.grid()

plt.savefig('myfig.png')