#Variable Declararion
L=18 #Latitude of earth station(degrees)
PE=-73 #Longitude of earth station(degrees)
Pss=-105 #Satellite location(degrees)
aGSO=42164 #Circumference of earth (km)
R=6371 #Radius of earth(km)
import math
#Calculation
B=PE-Pss #Angle between the planes containing a and c (degrees)
Rx=R*math.cos(L*3.142/180)*math.cos(B*3.142/180) #Geocentric-equitorial coordinate(km)
Ry=R*math.cos(L*3.142/180)*math.sin(B*3.142/180) #Geocentric-equitorial coordinate(km)
Rz=R*math.sin(L*3.142/180) #Geocentric-equitorial coordinate(km)
import numpy as np
r=np.array([Rx,Ry,Rz]) #Coordinates for local gravity direction
k=np.array([Rx-aGSO,Ry,Rz])#geocentric-equitorial coordinates for propagation direction
e=np.array([0,0,1]) # #geocentric-equitorial coordinates for polarization vector
f=np.cross(k,r) #Direction of normal to reference plane
modf=(f[0]**2+f[1]**2+f[2]**2)**0.5
g=np.cross(k,e)# Direction of normal to plane contaning e and k
h=np.cross(g,k) #Direction of polarization of the plane
modh=(h[0]**2+h[1]**2+h[2]**2)**0.5
p=(h/modh)
p[0]=round(p[0],3)
p[1]=round(p[1],3)
p[2]=round(p[2],3)
E=round(math.asin(np.dot(p,f)/modf)*180/3.142,2)
print "The Angle of polarization at given location is",E,"degrees"