%matplotlib inline
from IPython.display import Image

import numpy as np

"""
Explanation: 2014 final -- solutions
do detailed, labeled breakdown of layer question and radiance question
End of explanation
"""


np.searchsorted([2,4,5,6,10,12],[-2,-5,10,1,4,6,0,4,14], 'right')

"""
Explanation: Q1  -- searchsorted
a) (3)  Explain as clearly as you can what numpy's searchsorted function does.  Include your best
guess about what the variable vals contains after we execute statement [9].
Give an example of  how searchsorted is used in python programs we've encountered
in class -- you don't have to regurgitate the code, just describe why it is needed and what
specifically it accomplishes in one or more of our programs.
End of explanation
"""


from numpy import pi,exp,cos
import matplotlib.pyplot as plt


plt.close('all')
thetime=np.arange(0.,2.5*pi,0.05)
thewave=thetime
four5=45.*pi/180.
thirty=30*pi/180.
sixty=2.*thirty
onefifty=5.*thirty


fig1,axis1=plt.subplots(1,1)
#
#  pulse one -- phase shift is 45 degrees
#
axis1.plot(thewave,np.cos(thetime + four5),'k-',lw=5,label='pulse 1')
#
# pulse 2 -- phase shift is 150 degrees
#
axis1.plot(thewave,np.cos(thetime + onefifty),'k+',lw=5,label='pulse 2')
axis1.set_xlabel('horizontal position (one wavelength=2pi)')
axis1.set_ylabel('amplitude')
axis1.set_title(' ')
pos=[0.,0.25*pi,0.5*pi,0.75*pi,1.*pi,1.25*pi,1.5*pi,1.75*pi,2*pi,2.25*pi]
labels=['0','pi/4','pi/2','3pi/4','pi','5pi/4','6pi/4','7pi/4','2pi','9pi/4']
axis1.grid()
axis1.set_xticks(pos)
axis1.set_xticklabels(labels)
axis1.legend(loc='best')
fig1.savefig('phase_shift_final.png')

#pulse1
I=0.707  #cos phi is 0.707
Q= -(0.707) # =  sin phi is  0.707
#
#  so pulse 1 is in the upper right quadrant
#
#import math
#math.atan(0.707/0.707)*180./math.pi
#
# so upper right quadrant, 45 degrees
#
#pulse2
I=-0.866  #cos phi is -0.866
Q=(-0.5) #= sin phi is +0.5
#
# pulse 2 is in the upper left quadrant
#

#math.atan(0.5/0.866)*180./math.pi = 30 degrees
#first angle= +45
#pulse2=150 degrees
mrmax=0.1*600./4.
print("Maximum resolved velocity is {} m/s".format(mrmax))

#smallest angle is 150 - 45.
angle1=150 - 45.  #positive phase shift
vel1= -angle1/180.*mrmax  #-8.75 m/s into the radar
#
# smallest angle is positive (counterclockwise) so velocity is
# negative (into the radar)
# 
#
# for next angle, go clockwise
#
#second guess
angle2=360. - angle1
vel2=(angle2)/180.*mrmax #away from the radar
print(("first guess angle is {angle1:} degrees\n"
       "first guess velocity is {vel1:5.2f} m/s\n"
       "second guess angle is {angle2:} degrees\n"
       "second guess velocity is {vel2:5.2f} m/s").format_map(locals()))

"""
Explanation: Q2 - Derivation of Schwartzchild equation
See the Temperature changing with height section of the Schwartzchild notes
Q3 - Pulse Pair
A cloud radar operates at $\lambda=10\ cm$ with a PRF of 600 $s^{-1}$.
This figure shows the in-phase (coherence) plots of two wavetrains returning from a pulse pair
separated by 1/600th of a second.
End of explanation
"""


Image('figures/two_layer.png')

"""
Explanation: Q4
The figure below shows a nocturnal, two layer atmosphere over a black surface.  Given a mass absorption coefficient $\kappa = 0.001\ m^{2}/kg$
End of explanation
"""


sigma=5.67e-8
tau2=0.001*1500
tau1=0.001*1000.
tautot=tau1 + tau2
temp2=260.  #K
temp1=280.  #K
tempsfc=300.
t1=exp(-5./3.*tau1)
t2=exp(-5/3.*tau2)
eps1=(1.-t1)
eps2=(1.-t2)
Es=sigma*tempsfc**4.
E1=eps1*sigma*temp1**4.
E2=eps2*sigma*temp2**4.
#
# top of layer 2
#
Es12= Es*t1*t2
E12 = E1*t2
Etop = -E2 - Es12 - E12
#
# bot of layer 2
#
Es1=Es*t1
Ebot =  -Es1 + E2 - E1
Enet=Etop - Ebot
print(("optical depth of layer 1 is {tau1:5.2f}\n"
       "optical depth of layer 2 is {tau2:5.2f}\n"
       "combined optical depth  is {tautot:5.2f}\n\n"
       "Net downward fluxes at top of layer 2 ignoring signs:\n"
       "    Es*t1*t2: {Es12:5.2f} W/m^2\n"
       "    E1*t2: {E12:5.2f} W/m^2\n"
       "    E2:    {E2:5.2f} W/m^2\n"
       "    Etop = -E2 - Es12 - E12: {Etop:5.2f} W/m^2\n"
       "Net downward fluxes at base of layer 2:\n"
       "    Es*t1: {Es1:5.2f} W/m^2\n"
       "    E2: {E2:5.2f} W/m^2\n"
       "    E1: {E1:5.2f} W/m^2\n"
       "    Ebot = -Es1 + E1 - E1: {Ebot:5.2f}\n"
       "Net flux divergence across layer: Etop - Ebot: {Enet:5.2f}").format_map(locals()))

"""
Explanation: Now calculate the heating rate of layer 2 (the top layer).  Let upward fluxes be negative and downward fluxes be
positive
End of explanation
"""


Qr = Enet/1500/1004.*3600.*24.
print('heating rate in K/day: {:5.2f}'.format(Qr))

"""
Explanation: So the top layer is cooling at a rate of:
End of explanation
"""


B1=7  #W/m^2/micron/sr  -- 280 K
B2=4.9  #W/m^2/micron/sr  -- 260 K
angle=30.*pi/180.
#angle=60.*pi/180.
tau2=1.5/cos(angle)
tau1=1./cos(angle)
tr1=exp(-tau1)
tr2=exp(-tau2)
eps1=(1. - tr1)
eps2=(1. - tr2)

L=eps2*B2*tr1 + eps1*B1
print('radiance in W/m^2/micron/sr: {:5.2f}'.format(L))

del_omega=0.01
del_lambda=2
flux=L*del_omega*del_lambda
print('flux in W/m^2: {:5.2f}'.format(flux))

flux

"""
Explanation: Q5 radiance
For the figure in problem 4, suppose you were looking upward with an infrared instrument at an angle of 30 degrees off vertical.  Assuming that
the mass absorption coefficient is again $\kappa = 0.001\ m^{2}/kg$ for longwave photons, the telescope measures all wavelengths
between $10\ \mu m  - 12\ \mu$ and the telescope's field of view is 0.01 sr, find:
a) (4)  The radiance observed by the telescope
b) (4) The flux observed by the telescope
c) (2) The brightness temperature observed by the telescope
d) (2) Suppose you point the telescope closer to the horizon, at 60 degrees off zenith.  Does the brightness temperature
increase or decrease?  Explain.
End of explanation
"""


