{
 "metadata": {
  "name": "",
  "signature": "sha256:9c5b783b7af2e7759af05c065786293bfda0657b22d7595e2418aaa43901e7b0"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Chapter 2: Orbits and Launching Methods"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.1, Page 23"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "\n",
      "#Variable Declaration\n",
      "\n",
      "u=3.986*(10**14)         #Earth's Gravitational constant(m^3/sec^2)\n",
      "\n",
      "#Calculation\n",
      "\n",
      "n=(2*3.14)/(24*60*60)    #Mean Motion(rad/sec)\n",
      "a=((u/n**2)**(0.33333))/1000      #Radius of the orbit by kepler's 3rd law(km)\n",
      "a=round(a,2)\n",
      "\n",
      "#Result\n",
      "\n",
      "print \"The Radius of the circular orbit with 1 day period is:\", a,\"km\"\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The Radius of the circular orbit with 1 day period is: 42247.94 km\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.2, Page 29"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "\n",
      "#Variable Declaration\n",
      "NN=14.22296917        #Mean Motion (1/day)\n",
      "u=3.986005*(10**14)   #Earth's Gravitational COnstant(m^3/sec^2)\n",
      "\n",
      "#Calculation\n",
      "\n",
      "n0=(NN*2*3.142)/(24*60*60)                  #Mean Motion(rad/sec)\n",
      "a=((u/n0**2)**(0.33333))/1000   #Radius of the orbit by kepler's 3rd law(km)\n",
      "a=round(a,2)\n",
      "#Result\n",
      "\n",
      "print \"The Semimajor axis for given satellite parameters is:\", a,\"km\"\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The Semimajor axis for given satellite parameters is: 7193.97 km\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.3, Page 29"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "\n",
      "#Variable Declaration\n",
      "\n",
      "R=6371   #Mean Earth's radius(km)\n",
      "e=0.0011501  #Eccentricity\n",
      "a=7192.3     #Semimajor axis(km)\n",
      "\n",
      "#Calculation\n",
      "\n",
      "ra=a*(1+e)   #Radius Vector at apogee(km)\n",
      "rp=a*(1-e)   #Radius Vector at perigee(km)\n",
      "ha=round(ra-R,2)      #Apogee height(km)\n",
      "hp=round(rp-R,2)      #Perigee height(km)\n",
      "\n",
      "\n",
      "#Result\n",
      "\n",
      "print \"The Apogee height for given orbital parameters is:\", ha,\"km\"\n",
      "print \"The Apogee height for given orbital parameters is:\", hp,\"km\"\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The Apogee height for given orbital parameters is: 829.57 km\n",
        "The Apogee height for given orbital parameters is: 813.03 km\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.4, Page 31"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "\n",
      "#Variable Declaration\n",
      "\n",
      "aE=6378.141   #Earth's equitorial radius(km)\n",
      "e=0.002  #Eccentricity\n",
      "p=12     #period from perigee to perigee (hours)\n",
      "K1=66063.1704  #Constant (km^2)\n",
      "u=3.986005*(10**14)         #Earth's Gravitational constant(m^3/sec^2)\n",
      "\n",
      "\n",
      "#Calculation\n",
      "\n",
      "n=(2*math.pi)/(12*60*60)  #Mean Motion(rad/sec)\n",
      "anp=((u/n**2)**(0.3333))/1000      #Radius of the orbit by kepler's 3rd law(km)\n",
      "k2=(1-e**2)**1.5\n",
      "# Solving for perturbed value of semimajor axis\n",
      "import scipy\n",
      "import scipy.optimize\n",
      "def f(a):\n",
      "    y=(n-((u/a**3)**0.5)*(1+K1/a**2*k2))\n",
      "    return y\n",
      "a=scipy.optimize.fsolve(f,2)\n",
      "a=a/1000  #Converting a into km\n",
      "\n",
      "#Result\n",
      "\n",
      "print \"The nonperturbed value of semimajor axis is\", anp,\"km\"\n",
      "print \"The perturbed value of semimajor axis is\", a,\"km\"\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The nonperturbed value of semimajor axis is 26564.7679852 km\n",
        "The perturbed value of semimajor axis is [ 26610.22410209] km\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.5, Page 34"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "\n",
      "#Variable Declaration\n",
      "\n",
      "i=98.6328      #Angle(degrees)\n",
      "e=0.0011501    #eccentricity\n",
      "n=14.23304826  #Mean Motion(1/day)\n",
      "a=7192.3       #Semimajor axis(km)\n",
      "K1=66063.1704  #Known constant(km^2)\n",
      "\n",
      "#Calculation\n",
      "\n",
      "n0=(2*180*n)  #Mean Motion (deg/sec)\n",
      "K=(n0*K1)/((a**2)*((1-e**2)**2))  #Constant (deg/day)\n",
      "w=-K*math.cos(i*3.142/180)    #Rate of regression of nodes(deg/day)\n",
      "W=K*(2-2.5*(math.sin(i*3.142/180))**2)   #Rate of rotation of line of apsides(deg/day)\n",
      "w=round(w,3)\n",
      "W=round(W,3)\n",
      "#Results\n",
      "\n",
      "print \"The rate of regression of nodes is:\", w,\"deg/day\"\n",
      "print \"The rate of rotation of line of apsides is:\", W,\"deg/day\"\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The rate of regression of nodes is: 0.984 deg/day\n",
        "The rate of rotation of line of apsides is: -2.902 deg/day\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.6, Page 34"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Variable Declaration\n",
      "\n",
      "w=0.982  #rate of regression of nodes from Example 2.5(deg/day)\n",
      "W=-2.903 #rate of rotation of line of apsides from Example 2.5)deg/day)\n",
      "n=14.23304826  #Mean Motion(1/day)\n",
      "W0=113.5534   #Argument of perigee(deg)\n",
      "w0=251.5324   #Right ascension of the ascending node(deg)\n",
      "\n",
      "#Calculation\n",
      "\n",
      "PA=1/n   #Period \n",
      "w=round(w0+w*PA,2)  #New value of rate of regression of nodes(deg)\n",
      "W=round(W0+W*PA,2)  #New Value of rate of rotation of line of apsides(deg)\n",
      "\n",
      "#Result\n",
      "\n",
      "print \"New value of rate of regression of nodes is:\", w,\"deg\"\n",
      "print \"New value of rate of rotation of line of apsides is:\", W,\"deg\"\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "New value of rate of regression of nodes is: 251.6 deg\n",
        "New value of rate of rotation of line of apsides is: 113.35 deg\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.7, Page 38"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Calculation\n",
      "\n",
      "ndays=400*365    #Nominal number of days in 400years\n",
      "nleapyrs=400/4   #Nominal number of leap years\n",
      "gregoriandays=ndays+nleapyrs-3  #number of days in 400 years of Gregorian calendar\n",
      "gregavg=round(gregoriandays/float(400),2)    #number of days in 400 years of Gregorian calendar\n",
      "\n",
      "#Result\n",
      "print gregoriandays\n",
      "print \"The average length of the civil year in gregorian calender is:\", gregavg,\"days\"\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "146097\n",
        "The average length of the civil year in gregorian calender is: 365.24 days\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.8, Page 38"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Calculation and Results\n",
      "\n",
      "if 1987%4==0:\n",
      "  print \"1987 is a leap year\"\n",
      "else:\n",
      "  print \"1987 is not a leap year\"\n",
      "\n",
      "\n",
      "if 1988%4==0:\n",
      "  print \"1988 is a leap year\"\n",
      "else:\n",
      "  print \"1988 is not a leap year\"\n",
      "\n",
      "\n",
      "if 2000%400==0:\n",
      "  print \"2000 is a leap year\"\n",
      "else:\n",
      "  print \"2000 is not a leap year\"\n",
      "\n",
      "\n",
      "if 2100%400==0:\n",
      "  print \"2100 is a leap year\"\n",
      "else:\n",
      "  print \"2100 is not a leap year\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1987 is not a leap year\n",
        "1988 is a leap year\n",
        "2000 is a leap year\n",
        "2100 is not a leap year\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.9, Page 38"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Calculation\n",
      "days=324   #Number of days\n",
      "hours=math.floor(24*0.95616765)  # Number of hours\n",
      "decimalfraction1=24*0.95616765-hours\n",
      "minutes=math.floor(60*decimalfraction1)  # Number of minutes\n",
      "decimalfraction2=60*decimalfraction1-minutes\n",
      "seconds=round(60*decimalfraction2,2)  # Number of seconds\n",
      "\n",
      "#Result\n",
      "\n",
      "print decimalfraction1,decimalfraction2\n",
      "print \"An Epoch day has\", days,\"days\",hours,\"hours\",minutes,\"minutes\",seconds,\"seconds\"\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.9480236 0.881416\n",
        "An Epoch day has 324 days 22.0 hours 56.0 minutes 52.88 seconds\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.10, Page 39"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "\n",
      "#Variable Declaration\n",
      "\n",
      "y=2000      #year\n",
      "mon=12      #month\n",
      "dy=18       #day\n",
      "hours=13    #hours of the day\n",
      "minutes=0   #Minutes of the day\n",
      "seconds=0   #Seconds of the day\n",
      "\n",
      "\n",
      "#Calculation\n",
      "d=dy+(hours/24)+(minutes/(24*60))+seconds  #Days in December \n",
      "if mon<=2:\n",
      "   y=y-1\n",
      "   mon=mon+12\n",
      "else:\n",
      "   y=y\n",
      "   mon=mon\n",
      "\n",
      "A=math.floor(y/100)      #Converting years to days\n",
      "B=2-A+math.floor(A/4)    #Converting years to days\n",
      "C=math.floor(365.25*y)     #rounding the days \n",
      "D=math.floor(30.6001*(mon+1)) #Converting months to days\n",
      "JD=B+C+D+d+1720994.5     #Adding reeference to number of days\n",
      "\n",
      "\n",
      "#Result\n",
      "\n",
      "print \"The Julian day of given day is:\", JD,\"Days\"\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The Julian day of given day is: 2451896.5 Days\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.11, Page 41"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Variable Declaration\n",
      "\n",
      "JDref=2415020   #Reference Julian days\n",
      "JC=36525\n",
      "JD=2451897.0417 #Julian days with reference from Example 2.10\n",
      "\n",
      "#Calculation\n",
      "\n",
      "T=round((JD-JDref)/JC,4)  #Time in julian Centuries\n",
      "\n",
      "#Result\n",
      "\n",
      "print \"The time for given date is:\", T,\"Julian Centuries\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The time for given date is: 1.0096 Julian Centuries\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.12, Page 43"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Variable Declaration\n",
      "\n",
      "n=14.23304826  #Mean Motion (rev/day)\n",
      "M0=246.6853    #Mean Anomaly (degrees)\n",
      "t0=223.79688452 #Time of anomaly\n",
      "\n",
      "#Calculation\n",
      "\n",
      "T=round(t0-(M0/(n*360)),3)  #Time of perigee passage\n",
      "\n",
      "#Result\n",
      "\n",
      "\n",
      "print \"The time of perigee passage for NASA elements is:\", T,\"days\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The time of perigee passage for NASA elements is: 223.749 days\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.13, Page 44"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#Variable Declaration\n",
      "\n",
      "M=205  #Mean anomaly(degrees)\n",
      "e=0.0025  #Eccentricity\n",
      "E=math.pi #Initial guess value for eccentric anomaly\n",
      "\n",
      "#Calculation\n",
      "\n",
      "import scipy\n",
      "import scipy.optimize\n",
      "def f(E):\n",
      "    y=M-E+e*math.sin(E)\n",
      "    return y\n",
      "E=scipy.optimize.fsolve(f,3.142)\n",
      "E=round(E,3)\n",
      "\n",
      "print \"The Eccentric anomaly is:\",E,\"degrees\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The Eccentric anomaly is: 204.998 degrees\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.14, Page 45"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#Variable Declaration\n",
      "\n",
      "n=14.2171401*2*math.pi/86400  #Mean motion (rad/sec)\n",
      "M=204.9779+0.001*180*5/math.pi #Mean anomaly(rad)\n",
      "e=9.5981*10**-3   #Eccentricity\n",
      "a=7194.9   #Semimajor axis(km)\n",
      "\n",
      "#Calculation\n",
      "\n",
      "v=(M*math.pi/180)+2*e*math.sin(M*math.pi/180)+(5*e**2*math.sin(2*M*math.pi)/(4*180)) #True Anomaly (radians)\n",
      "v=v*180/math.pi  #True anomaly(degrees)\n",
      "r=a*(1-e**2)/(1+e*math.cos(v))  #Magnitude of radius vector after 5s(km)\n",
      "v=round(v,2)\n",
      "r=round(r,2)\n",
      "\n",
      "#Results\n",
      "\n",
      "print \"The true anomaly is:\",v,\"degrees\"\n",
      "print \"The magnitude of radius vector 5s after epoch is:\",r,\"km\"\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The true anomaly is: 204.79 degrees\n",
        "The magnitude of radius vector 5s after epoch is: 7252.02 km\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.15, Page 46"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#Variable Declaration\n",
      "\n",
      "v=204.81   #True anomaly(degrees) from Example 2.14\n",
      "r=7257     #Magnitude of radius vector(km) from Example 2.14\n",
      "\n",
      "#Calculation\n",
      "\n",
      "rP=r*math.cos(v*math.pi/180)  #P coordinate of radius vector(km)\n",
      "rQ=r*math.sin(v*math.pi/180)  #Q coordinate of radius vector(km)\n",
      "rP=round(rP,2)\n",
      "rQ=round(rQ,2)\n",
      "#Result\n",
      "\n",
      "print \"r in the perifocal coordinate system is\", rP,\"Pkm\", rQ,\"Qkm\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "r in the perifocal coordinate system is -6587.21 Pkm -3045.11 Qkm\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Exanple 2.17, Page 50"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#Variable Declaration\n",
      "\n",
      "T=1.009638376  #Time in Julian centuries from Example 2.11\n",
      "UT=13  #Universal time(hours)\n",
      "\n",
      "#Calculation\n",
      "\n",
      "GST=(99.6910+36000.7689*T+0.004*T**2)*3.142/180  #GST(radians)\n",
      "UT=2*math.pi*UT/24  #Universal time converted to fraction of earth rotation (radians)\n",
      "\n",
      "GST=GST+UT \n",
      "\n",
      "\n",
      "GST=(math.fmod(GST,2*math.pi))*180/math.pi#using fmod multiplr revolutions are removed (degrees)\n",
      "GST=round(GST,2)\n",
      "\n",
      "#Result\n",
      "\n",
      "print \"The GST for given date and time is \", GST,\"degrees\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The GST for given date and time is  287.18 degrees\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.18, Page 51"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#Variable Declararion\n",
      "\n",
      "WL=-89.26  #Expressing the longitude in degrees west\n",
      "GST=282.449  #GST from Example 2.17 (degrees)\n",
      "\n",
      "#Calculation\n",
      "\n",
      "EL=2*math.pi+WL  #Longitude in degrees East \n",
      "LST=(GST+EL)*math.pi/180       #LST(radians)\n",
      "LST=(math.fmod(LST,2*math.pi))*180/math.pi   #fmod removes multiple revolutions(Degrees)\n",
      "LST=round(LST,2)\n",
      "\n",
      "#Results\n",
      "\n",
      "print \"LST for Thunder Bay on given day is:\" ,LST,\"Degrees\"\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "LST for Thunder Bay on given day is: 199.47 Degrees\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.19, Page 52"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#Variable Declaration\n",
      "\n",
      "LST=167.475 #LST(degrees)\n",
      "LE=48.42    #Latitude at thunder bay(degrees)\n",
      "H=200       #Height above sea level(metres)\n",
      "aE=6378.1414 #Semimajor axis(km)\n",
      "eE=0.08182   #Eccentricity\n",
      "\n",
      "#Calculation\n",
      "\n",
      "l=(aE/math.sqrt(1-eE**2*math.sin(LE*3.142/180)**2)+H/1000)*math.cos(LE*3.142/180)  \n",
      "\n",
      "z=((aE*(1-eE**2))/math.sqrt(1-eE**2*math.sin(LE*3.142/180)**2)+H/1000)*math.sin(LE*3.142/180)  \n",
      "\n",
      "RI=round(l*math.cos(LST*3.142/180),2) #I component of radius vector at thunder bay(km)\n",
      "\n",
      "RJ=round(l*math.sin(LST*3.142/180),2) #J component of radius vector at thunder bay(km)\n",
      "\n",
      "RK=round(z,2)      #Z component of radius vector at thunder bay(km)\n",
      "\n",
      "R=math.sqrt(RI**2+RJ**2+RK**2)\n",
      "R=round(R,2)\n",
      "\n",
      "\n",
      "#Results\n",
      "\n",
      "print \"The Radius vector components are\" ,RI,\"ikm+\",RJ,\"jkm+\",RK,\"kkm\"\n",
      "print \"The Magnitude of radius component is\" ,R,\"km\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The Radius vector components are -4139.81 ikm+ 918.02 jkm+ 4748.46 kkm\n",
        "The Magnitude of radius component is 6366.21 km\n"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.20, Page 55"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Variable Declaration\n",
      "\n",
      "PI=-1280  #I component of range vector for a satellite(km)\n",
      "PJ=-1278  #J component of range vector for a satellite(km)\n",
      "PK=66     #K component of range vector for a satellite(km)\n",
      "GST=240   #GST(degrees)\n",
      "LE=48.42  #Latitude(Degrees)\n",
      "PE=-89.26 #Longitude(Degrees)\n",
      "H=200     #Height above mean sea level(metres)\n",
      "aE=6378.1414 #Semimajor axis(km)\n",
      "eE=0.08182   #Eccentricity\n",
      "\n",
      "import math\n",
      "#Calculation\n",
      "\n",
      "l=(aE/math.sqrt(1-eE**2*math.sin(LE*3.142/180)**2)+H/1000)*math.cos(LE*3.142/180)  \n",
      "\n",
      "z=((aE*(1-eE**2))/math.sqrt(1-eE**2*math.sin(LE*3.142/180)**2)+H/1000)*math.sin(LE*3.142/180)  \n",
      "\n",
      "SE=(math.atan(z/l))*180/3.142 #Geocentric latitude angle (degrees)\n",
      "LST=240+PE\n",
      "\n",
      "from numpy import matrix\n",
      "from numpy import linalg\n",
      "a=math.sin(SE*3.142/180)*math.cos(LST*3.142/180)\n",
      "b=math.sin(SE*3.142/180)*math.sin(LST*3.142/180)\n",
      "c=-math.cos(SE*3.142/180)\n",
      "d=-math.sin(LST*3.142/180)\n",
      "e=math.cos(LST*3.142/180)\n",
      "f=0\n",
      "g=math.cos(SE*3.142/180)*math.cos(LST*3.142/180)\n",
      "h=math.cos(SE*3.142/180)*math.sin(LST*3.142/180)\n",
      "i=math.sin(SE*3.142/180)\n",
      "\n",
      "D = matrix( [[a,b,c],[d,e,f],[g,h,i]]) \n",
      "\n",
      "P=matrix( [[PI],[PJ],[PK]] )\n",
      "\n",
      "R=D*P  #Components of range of earth station\n",
      "Ro=math.sqrt(R[0,0]**2+R[1,0]**2+R[2,0]**2) #Magnitude of range of earth station(km)\n",
      "Ro=round(Ro,2)  \n",
      "El=math.asin(R[2,0]/Ro) #Antenna elevation angle for the earth station(radians) \n",
      "El=round(El*180/3.142,2) #Converting El to degrees\n",
      "alpha=(math.atan(R[1,0]/R[2,0]))*180/3.142\n",
      "if (R[0,0]<0 and R[1,0]>0):\n",
      "   Aza=alpha\n",
      "else:\n",
      "   Aza=0\n",
      "if (R[0,0]>0 and R[1,0]>0):\n",
      "   Azb=180-alpha\n",
      "else:\n",
      "   Azb=0\n",
      "if (R[0,0]>0 and R[1,0]<0):\n",
      "   Azc=180+alpha\n",
      "else:\n",
      "   Azc=0\n",
      "if (R[0,0]<0 and R[1,0]<0):\n",
      "   Azd=360-alpha\n",
      "else:\n",
      "   Azd=0\n",
      "Az=round(Aza+Azb+Azc+Azd,2)  #Azimuth angle (degrees)\n",
      "\n",
      "print \"The magnitude of range of earth station is\" ,Ro,\"km\"\n",
      "print \"The antenna elevation angle for the earth station are\",El,\"degrees\"\n",
      "print \"The Azimuth angle for the earth station is\",Az,\"degrees\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The magnitude of range of earth station is 1809.98 km\n",
        "The antenna elevation angle for the earth station are 12.03 degrees\n",
        "The Azimuth angle for the earth station is 102.24 degrees\n"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 2.21, Page 58"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Variable Declaration\n",
      "\n",
      "rI=-4685.3  #I component of radius vector from Example 2.16(km)\n",
      "rJ=5047.7   #J component of radius vector from Example 2.16(km)\n",
      "rK=-3289.1  #K component of radius vector from Example 2.16(km)\n",
      "aE=6378.1414 #Semimajor axis (km)\n",
      "eE=0.08182   #Eccentricity\n",
      "\n",
      "#Calculation\n",
      "\n",
      "r=math.sqrt(rI**2+rJ**2+rK**2)\n",
      "a=math.pi  #Guess value for LST(radians)\n",
      "b=math.atan(rK/rI)   #Guess Value for latitude(radians)\n",
      "c=r-aE   #Guess value for height(km)\n",
      "\n",
      "from scipy.optimize import fsolve\n",
      "\n",
      "\n",
      "def equations(p):\n",
      "    L,h,LST = p\n",
      "    return (rI-((aE/math.sqrt(1-eE**2*math.sin(L)**2))+h)*math.cos(L)*math.cos(LST),rJ-((aE/math.sqrt(1-eE**2*math.sin(L)**2))+h)*math.cos(L)*math.sin(LST),rK-((aE*(1-eE**2)/math.sqrt(1-eE**2*math.sin(L)**2))+h)*math.sin(L))\n",
      "\n",
      "L,h,LST =  fsolve(equations, (b,c,a))\n",
      "L=round(L*180/3.142,2)  #Converting L into degrees\n",
      "h=round(h)\n",
      "LST=round(LST*180/3.142,1) #Converting LST into degrees\n",
      "\n",
      "print \"The latitude of subsatellite is\",L,\"degrees\"\n",
      "print \"The height of subsatellite is\",h,\"km\"\n",
      "print \"The LST of subsatellite is\",LST,\"degrees\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The latitude of subsatellite is -25.65 degrees\n",
        "The height of subsatellite is 1258.0 km\n",
        "The LST of subsatellite is 132.9 degrees\n"
       ]
      }
     ],
     "prompt_number": 23
    }
   ],
   "metadata": {}
  }
 ]
}