___  ____  ____  ____  ____(R)
                                                      /__    /   ____/   /   ____/   
                                                     ___/   /   /___/   /   /___/    
                                                       Statistics/Data Analysis      
      
      -------------------------------------------------------------------------------
            name:  <unnamed>
             log:  /Users/robertyaffee/Documents/data/research/chwk/phase3/geodesicDi
      > stanceComputations.smcl
        log type:  smcl
       opened on:  20 Jan 2012, 18:05:09
      
     1 . *  Robert Yaffee  20 Jan 2012    Haversine Geodesic distance computations
     2 . *  Chornobyl sequelae study
     3 . 
     4 . set more off
      
     5 . cap label var lat1 "Chornobyl latitude in decimal degrees"
      
     6 . cap label var lon1 "Chornobyl longitud in decimal degrees"
      
     7 . cap label var lat2 "1986 residential latitude"
      
     8 . cap label var lon2 "1986 residential longitude"
      
     9 . cap label var sinlat1 "sin of Chornobyl latitude in decimal degrees"
      
    10 . cap label var coslat1 "cos of Chornobyl latitude in decimal degrees" 
      
    11 . cap label var sinlat2 "sin of residential latitude in decimal degrees at time
      >  of accident"
      
    12 . cap label var coslat2 "cos of residential latitude in decimal degrees at time
      >  of accident"
      
    13 . cap label var cosdlon "cos of difference between latitudes of Ch and April 26
      >  1986 residence"
      
    14 . cap label var cosdlonrad "cosine of difference between latitudes of Ch and Ap
      > rl 26 1986 residence in radians"
      
    15 . cap label var sindlatr "sine of difference in Ch and April 26 1986 residentia
      > l latitude in radians"
      
    16 . cap label var sindlonr "sine of difference in Ch and April 26 1986 residentia
      > l longitude in radians"
      
    17 . 
    18 .  //  the objectiveis to compute the distance between Chornobyl
    19 .  //  and residences of respondents in 1986.  To do so,
    20 .  //  we will use the Haversine formula that has been applied
    21 .  //  in spherical trigonometry for this purpose.
    22 .  //    algorithmic summary
    23 .  //   1.  convert degrees and minutes to decimal degrees
    24 .  //   2.  convert decimal degrees to radians.
    25 .  //   3.  Compute the sins and cosine needed to find the distance
    26 .  //   4.  the distance will be found and must be reconverted to the original 
      > metric
    27 . 
    28 . 
    29 . 
    30 . *=================== Baseline converstion factors============================
      > ========
    31 . 
    32 . *------------------ 1 mile = 1.609344 km
    33 . *------------------ 1 km = .62137119 miles
    34 . *------------------ earth radius = 6371 km
    35 . *-------------------earth radius = 3963.1 miles
    36 . 
    37 . *------------------ decimal degrees = degrees + minutes/60    // seconds are 
      > ignored
    38 . *------------------ angle in radians = .0174329 * degrees    
    39 . 
    40 . 
    41 . *-------------------Chornobyl coordinates    latitude in degrees and minutes 
      > with 0 seconds =  51¡23"N   longitude= 30¡6" E
    42 . cap gen LatChrnDD = 51 + 23/60
      
    43 . cap gen LonChrndD = 30 + 6/60
      
    44 . cap gen lon1dd = 30.08
      
    45 . 
    46 . 
    47 . *------------------Kiev City Center Coordinates latitude in degree  = 50 27" 
      > N    longitude = 30 32" E
    48 . 
    49 . *------ Nonspherical map distance according to Google Maps 19 Jan 2012
    50 . * Baseline computation by google earth of distance from center of Kiev to Cho
      > rnobyl
    51 . * distance according to Google maps from Chornobyl to Kiev City 68.8 miles = 
      > 110.53 km
    52 . 
    53 . 
    54 . *============================================================================
      > ========================
    55 . 
    56 . 
    57 . //       Haversine formula  for great circle distances between geodesic coord
      > inates of lat and lon  
    58 . //    Sources: http://en.wikipedia.org/wiki/Great-circle_distance#Radius_for_
      > spherical_Earth,  accessed 18 20 Jan 2012
    59 . //     http://mathworld.wolfram.com/GreatCircle.html,  accessed 20 Jan 2012
    60 . //      http://maps.google.com
    61 . //      http://www.movable-type.co.uk/scripts/latlong-vincenty.html  by Chris
      >  Veness accessed 20 Jan 2012
    62 . //      Stata ado file written by Austin Nichols to compute Geodesic distance
      > s by the Thaddeus Vincenty formula, available
    63 . //         from Stata Software Components archive managed by Chris Baum at Bo
      > ston College. accessed 20 Jan 2012
    64 . 
    65 . *============================================================================
      > =======================
    66 . 
    67 .            
    68 . //      Distance = radius of earth( Radearth ) * Central angle betwen two poi
      > nts (CAR) measured in radians.
    69 . //  This distance is measured in km.
    70 . 
    71 . //      where
    72 . 
    73 . //  CA = acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(dlon) )      // C
      > entral angle distance with angle is measured in radians 
    74 . // Radians = angle in degrees * .01745329
    75 . // Distance =  CA * Radearth = 6371.01
    76 . 
    77 . 
    78 . replace RadEarthkm = 6371.01                       // in km 
      (0 real changes made)
      
    79 . label var RadEarthkm "Radius of the earth in km"
      
    80 . list RadEarthkm if _n < 5
      
           +----------+
           | RadEa~km |
           |----------|
        1. |  6371.01 |
        2. |  6371.01 |
        3. |  6371.01 |
        4. |  6371.01 |
           +----------+
      
    81 . replace RadEarthMi =  3963.1
      (0 real changes made)
      
    82 . label var RadEarthMi "Radius of earth in miles"
      
    83 . list RadEarthMi if _n < 5
      
           +----------+
           | RadEar~i |
           |----------|
        1. |   3963.1 |
        2. |   3963.1 |
        3. |   3963.1 |
        4. |   3963.1 |
           +----------+
      
    84 . 
    85 . replace lat2dd = latdacc + latmacc/60
      (0 real changes made)
      
    86 . replace lon2dd = londacc + lonmacc/60
      (0 real changes made)
      
    87 . list lat2dd if _n < 5
      
           +----------+
           |   lat2dd |
           |----------|
        1. | 50.73333 |
        2. | 50.73333 |
        3. | 50.73333 |
        4. | 50.73333 |
           +----------+
      
    88 . list lon2dd if _n < 5
      
           +----------+
           |   lon2dd |
           |----------|
        1. | 30.86667 |
        2. | 30.86667 |
        3. | 30.86667 |
        4. | 30.86667 |
           +----------+
      
    89 . 
    90 . // use angles in decimal degrees to begin computations  
    91 . cap gen lat1dd = LatChrnDD                      // Chornobyl latitude in deci
      > mal degrees
      
    92 . label var lat1dd "Latitude of Chorn. in decimal degrees"
      
    93 . cap gen lon1dd = LonChrnDD                      // Chornobyl longitude in dec
      > imal degrees
      
    94 . label var lon1dd "Longitude of Chornobyl in decimal degrees"
      
    95 . cap rename Latdd1986 Lat1986DD                  //  latitude location at time
      >  of accident on part of respondent
      
    96 . cap rename Londd1986 Lon1986DD                  //  longitude location at tim
      > e of accident  "
      
    97 . cap gen lat2dd = Lat1986DD                     //  latitude   "
      
    98 . cap gen lon2dd = Lon1986DD                     //  longitude  "
      
    99 . label var lat2dd "Latitude of respondent  in decimal degrees in April 1986"
      
   100 . label var lon2dd "Latitude of respondent in decimal degrees in April 1986"
      
   101 . // we convert these angle to radians for use in trig formulae
   102 . cap gen  lat1r =LatChrnDD*.01745329                       // angle in radians
      
   103 . cap gen  lon1r =LonChrnDD*.01745329                       // angle in radians
      
   104 . cap gen  lat2r =lat2dd*.01745329                          // angle in radians
      
   105 . cap gen  lon2r =Lon2dd*.01745329                          // angle in radians
      
   106 . 
   107 . // we substitute radians for degrees in trig formulae
   108 . cap gen  sinlat1r = sin(lat1r)                     //   sin of Chornobyl lati
      > tude in radians
      
   109 . cap gen  coslon1r = cos(lon1r)                     //   cos of Chornobyl long
      > itude in radians
      
   110 . cap gen  sinlat2r= sin(lat2r)                      // sin of 1986 latitude
      
   111 . cap gen  coslat2r = cos(lat2r)                     // cos of 1986 latitude
      
   112 . cap gen  cosdlonr = cos(lon2r-Lon1r)       // cos of longitudinal difference 
      > between the two in radians
      
   113 . cap gen sindlatr = sin(lat2r-lat1r)        // sin of latitudinal difference b
      > etween the two in radians
      
   114 . cap gen sindlonr = sin(lon2r - lon1r)
      
   115 . cap gen dlatr = (lat1dd - lat2dd)*01745329   
      
   116 . label var dlatr "Difference in latitudes measured in radians between Ch and 1
      > 986 april res location"
      note: label truncated to 80 characters
      
   117 . 
   118 . *---------------------  2 formulae for central angle ------------------------
      > -------------------------------
   119 . 
   120 . replace dx21 = 2 * asin((sin(dlatr * .5))^2 + cos(lat1r) * cos(lat2r)*(sin(dl
      > onr*.5))^2)       // formula 1
      (0 real changes made)
      
   121 . replace DX21 = acos(sin(lat1r)*sin(lat2r) + cos(lat1r)*cos(lat2r)*cos(dlonr) 
      > )                 // formula 2
      (0 real changes made)
      
   122 . label var dx21 "Dist. from point 1 to point 2 in km divided by the earth radi
      > us in km" 
      
   123 . 
   124 . 
   125 . replace DX_Chkm = DX21 * 6371.01   // distance in km
      (0 real changes made)
      
   126 . replace DX_Chmi = DX21 * 3963.1      // distance in miles
      (0 real changes made)
      
   127 . 
   128 . 
   129 . replace dxfmChkm = dx21 * 6371.01   // distance in km
      (698 real changes made)
      
   130 . replace dxfmChMi = dx21 * 3963.1    // distance in miles
      (698 real changes made)
      
   131 . label var dxfmChkm "Distance from Chornobyl in km via Haversine formula 1"
      
   132 . 
   133 . 
   134 . 
   135 . 
   136 . list DX_Chkm dxfmChkm DX_Chmi dxfmChMi sett1r1 if _n < 10
      
           +-----------------------------------------------------+
           |  DX_Chkm   dxfmChkm    DX_Chmi   dxfmChMi   sett1r1 |
           |-----------------------------------------------------|
        1. | 74.90189   19457.48   46.59288   12103.57      Kyiv |
        2. | 74.90189   19457.48   46.59288   12103.57      Kyiv |
        3. | 74.90189   19457.48   46.59288   12103.57      Kyiv |
        4. | 74.90189   19457.48   46.59288   12103.57      Kyiv |
        5. | 74.90189   19457.48   46.59288   12103.57      Kyiv |
           |-----------------------------------------------------|
        6. | 80.79928   19462.79   50.26136   12106.87      Kyiv |
        7. | 22.50504   1851.379    13.9993   1151.654   Pripyat |
        8. | 74.90189   19457.48   46.59288   12103.57      Kyiv |
        9. | 74.90189   582.7817   46.59288   362.5205      Kyiv |
           +-----------------------------------------------------+
      
   137 . 
   138 . * Method 2   Computation of great circle distance from Chornobyl via Haversin
      > e formula
   139 . 
   140 . 
   141 . //  Haversine formula uses 
   142 . //    1. dms converted to decimal degrees            lataccDD lonaccDD
   143 . //    2. decimal degrees converted to radians        lat2r lon2r
   144 . //    3. differences in latitude longitude measured in radians  dlatr and dlo
      > nr
   145 . 
   146 . cap gen lataccDD = latdacc + (latmacc/60)
      
   147 . label var lataccDD "latitude at time of Chorn. accident in decimal degrees"
      
   148 . cap gen lonaccDD = londacc + (lonmacc/60)
      
   149 . label var lonaccDD "longitude at time of Chorn. accident in decimal degrees"
      
   150 . cap gen lat2r =LataccDD*.01745329                       // angle in radians
      
   151 . cap gen lon2r =LonaccDD*.01745329                       // angle in radians
      
   152 . label var lat2r "Latitude in radians at time of Chorn. accident"
      
   153 . label var lon2r "Longitude in radians at time of Chorn. accident"
      
   154 . 
   155 . 
   156 . cap gen dlatr = lat2r - lat1r 
      
   157 . cap gen dlonr = lon2r - lon1r
      
   158 . label var dlatr "Difference in latitude in radians" 
      
   159 . label var dlonr "Difference in longitude in radians"
      
   160 . 
   161 . label var sinlat1r "sin of latitude of Chornobyl in radians"
      
   162 . label var coslat1r "cos of latitutde of Chornobyl in radians"
      
   163 . label var coslon1r "cos of longitude of Chornobyl in radians"
      
   164 . label var dlatr "Latitudinal difference measured in radians  bet Chrn. & resi
      > dence"
      
   165 . 
   166 . label var dlonr "Longitudinal difference betw. Chrn. & residence in radians"
      
   167 . label var dlonr "Longitudinal difference betw. Chrn. & residence in radians"
      
   168 . label var sinlat2r "Sin of latitude of residence in radians at time of Chrn. 
      > accident"
      
   169 . label var coslat2r "Cos of latitude of residence in radians at time of Chrn. 
      > accident"
      
   170 . 
   171 . 
   172 . //  arc span 
   173 . replace Thearc = ((sin(dlatr/2))^2 + cos(lat1r)*cos(lat2r)*(sin(dlonr/2))^2)
      (0 real changes made)
      
   174 . label var Thearc "Harversine arc computation part 1"
      
   175 . 
   176 . //given
   177 . 
   178 . *----------------------------------------------------------------------// for
      > mula 3 for central angle
   179 . replace ca3 = 2*atan2(sqrt(Thearc), sqrt(1-Thearc))
      (0 real changes made)
      
   180 . label var ca3 "Central angle computation for Haversine formula in radians"
      
   181 . 
   182 . 
   183 . replace dxfmChkm = 6371.01 * ca3
      (698 real changes made)
      
   184 . label var dxfmChkm "Residential distance from Chornobyl in km"
      
   185 . 
   186 . replace RadEarthMi = 3963.1
      (0 real changes made)
      
   187 . label var RadEarthMi "Radius of earth in miles"
      
   188 . 
   189 . replace dxfmChMi = 3963.1* ca3
      (698 real changes made)
      
   190 . label var dxfmChMi "Residential dist from Chornobyl in miles"
      
   191 . 
   192 . replace HavKm = dxfmChkm    // Haversine formula used to compute distance fro
      > m Chornobyl in km
      (0 real changes made)
      
   193 . replace Havmil = dxfmChmi   // Haversine formula used to compute distance fro
      > m Chornobyl in miles
      (0 real changes made)
      
   194 . label var HavKm "Distance from Chornobyl in Km"
      
   195 . label var Havmil "Distance from Chornobyl in miles"
      
   196 . list sett1r1 HavKm Havmil DX_Chkm dxfmChkm DX_Chmi dxfmChMi sett1r1 if _n < 1
      > 0
      
           +-------------------------------------------------------------------------
      ------------+
           | sett1r1      HavKm     Havmil    DX_Chkm   dxfmChkm    DX_Chmi   dxfmChM
      > i   sett1r1 |
           |-------------------------------------------------------------------------
      ------------|
        1. |    Kyiv   19620.78   46.59288   74.90189   19620.78   46.59288   12205.1
      > 4      Kyiv |
        2. |    Kyiv   19620.78   46.59288   74.90189   19620.78   46.59288   12205.1
      > 4      Kyiv |
        3. |    Kyiv   19620.78   46.59288   74.90189   19620.78   46.59288   12205.1
      > 4      Kyiv |
        4. |    Kyiv   19620.78   46.59288   74.90189   19620.78   46.59288   12205.1
      > 4      Kyiv |
        5. |    Kyiv   19620.78   46.59288   74.90189   19620.78   46.59288   12205.1
      > 4      Kyiv |
           |-------------------------------------------------------------------------
      ------------|
        6. |    Kyiv   19624.54   50.26136   80.79928   19624.54   50.26136   12207.4
      > 8      Kyiv |
        7. | Pripyat   4973.792    13.9993   22.50504   4973.792    13.9993   3093.95
      > 8   Pripyat |
        8. |    Kyiv   19620.78   46.59288   74.90189   19620.78   46.59288   12205.1
      > 4      Kyiv |
        9. |    Kyiv   2745.761   46.59288   74.90189   2745.761   46.59288   1708.00
      > 6      Kyiv |
           +-------------------------------------------------------------------------
      ------------+
      
   197 . 
   198 . *********************  Computation of great circle dx from Chornobyl at time 
      > of accident via Spherical law of cosines
   199 . 
   200 . //  Spherical law of cosines    
   201 . 
   202 . 
   203 . cap drop SLOCDX
      
   204 . cap gen SLOCDXkm =acos(sin(lat1r)*sin(lat2r)+cos(lat1r)*cos(lat2r)*cos(dlonr)
      > )*6371.01   
      
   205 . label var SLOCDXkm "Spherical law of cosines computation of dx from Chornobyl
      >  in km"
      
   206 . label var SLOCDXMi "Spherical law of cosines computation of dx from Chornobyl
      >  in Miles"
      
   207 . 
   208 . 
   209 . list sett1r1 HavKm Havmil SLOCDXkm if _n < 10 
      
           +------------------------------------------+
           | sett1r1      HavKm     Havmil   SLOCDXkm |
           |------------------------------------------|
        1. |    Kyiv   19620.78   46.59288   74.90189 |
        2. |    Kyiv   19620.78   46.59288   74.90189 |
        3. |    Kyiv   19620.78   46.59288   74.90189 |
        4. |    Kyiv   19620.78   46.59288   74.90189 |
        5. |    Kyiv   19620.78   46.59288   74.90189 |
           |------------------------------------------|
        6. |    Kyiv   19624.54   50.26136   80.79928 |
        7. | Pripyat   4973.792    13.9993   22.50504 |
        8. |    Kyiv   19620.78   46.59288   74.90189 |
        9. |    Kyiv   2745.761   46.59288   74.90189 |
           +------------------------------------------+
      
   210 . 
   211 . 
   212 . 
   213 . cap gen SLOCDXMi = SLOCDXkm/1.609
      
   214 . 
   215 . list id sett1r1 lat1dd lon1dd latdacc latmacc londacc lonmacc HavKm Havmil if
      >  sett1r1=="Pripyat"
      
           +-------------------------------------------------------------------------
      ----------------------+
           |  id   sett1r1   lat1dd   lon1dd   latdacc   latmacc   londacc   lonmacc 
      >      HavKm     Havmil |
           |-------------------------------------------------------------------------
      ----------------------|
        7. |  22   Pripyat     51.4    30.08        51        10        30        30 
      >   4973.792    13.9993 |
       40. | 167   Pripyat     51.4    30.08        51        49        30        11 
      >   844.4633   37.70615 |
      324. | 767   Pripyat     51.4    30.08        51        40        30         8 
      >   276.8539   27.56934 |
      375. |  64   Pripyat     51.4    30.08        51        10        30        30 
      >   4973.792    13.9993 |
           +-------------------------------------------------------------------------
      ----------------------+
      
   216 . 
   217 . 
   218 . list id sett1r1 lat1dd lon1dd latdacc latmacc londacc lonmacc HavKm Havmil  i
      > f sett1r1=="Ivankiv"
      
           +-------------------------------------------------------------------------
      ----------------------+
           |  id   sett1r1   lat1dd   lon1dd   latdacc   latmacc   londacc   lonmacc 
      >      HavKm     Havmil |
           |-------------------------------------------------------------------------
      ----------------------|
      143. | 542   Ivankiv     51.4    30.08        50        31        31         6 
      >   6925.673   64.73348 |
      583. | 438   Ivankiv     51.4    30.08        50        31        31         6 
      >   6925.673   64.73348 |
           +-------------------------------------------------------------------------
      ----------------------+
      
   219 . list id sett1r1 dxfmChkm dxfmChmi SLOCDXkm SLOCDXMi HavKm Havmi if _n < 20
      
           +-------------------------------------------------------------------------
      -----------+
           | id      sett1r1   dxfmChkm    dxfmChmi   SLOCDXkm   SLOCDXMi      HavKm 
      >     Havmil |
           |-------------------------------------------------------------------------
      -----------|
        1. |  3         Kyiv   19620.78   46.592876   74.90189   46.55183   19620.78 
      >   46.59288 |
        2. | 13         Kyiv   19620.78   46.592876   74.90189   46.55183   19620.78 
      >   46.59288 |
        3. | 14         Kyiv   19620.78   46.592876   74.90189   46.55183   19620.78 
      >   46.59288 |
        4. | 15         Kyiv   19620.78   46.592876   74.90189   46.55183   19620.78 
      >   46.59288 |
        5. | 16         Kyiv   19620.78   46.592876   74.90189   46.55183   19620.78 
      >   46.59288 |
           |-------------------------------------------------------------------------
      -----------|
        6. | 20         Kyiv   19624.54   50.261356   80.79928   50.21708   19624.54 
      >   50.26136 |
        7. | 22      Pripyat   4973.792   13.999305   22.50504   13.98697   4973.792 
      >    13.9993 |
        8. | 36         Kyiv   19620.78   46.592876   74.90189   46.55183   19620.78 
      >   46.59288 |
        9. | 40         Kyiv   2745.761   46.592876   74.90189   46.55183   2745.761 
      >   46.59288 |
       10. | 56         Kyiv   19620.78   46.592876   74.90189   46.55183   19620.78 
      >   46.59288 |
           |-------------------------------------------------------------------------
      -----------|
       11. | 58         Kyiv   11018.24   46.592876   74.90189   46.55183   11018.24 
      >   46.59288 |
       12. | 62    Chornobyl   18965.11   .54694138   .8792534   .5464596   18965.11 
      >   .5469414 |
       13. | 69     Vishneve    2468.62   57.205126   91.96195   57.15472    2468.62 
      >   57.20512 |
       14. | 71      Boyarka   10741.33   52.851079   84.96247   52.80452   10741.33 
      >   52.85108 |
       15. | 75      Boyarka   16319.37   53.666692   86.27364   53.61941   16319.37 
      >   53.66669 |
           |-------------------------------------------------------------------------
      -----------|
       16. | 78     Polesske   3589.894   13.999305   22.50504   13.98697   3589.894 
      >    13.9993 |
       17. | 82   Barishevka   2470.347    80.57438     129.53   80.50339   2470.347 
      >   80.57438 |
       18. | 84         Kyiv   19620.78   46.592876   74.90189   46.55183   19620.78 
      >   46.59288 |
       19. | 85         Kyiv   19620.78   46.592876   74.90189   46.55183   19620.78 
      >   46.59288 |
           +-------------------------------------------------------------------------
      -----------+
      
   220 . save chwide19jan2012, replace
      file chwide19jan2012.dta saved
      
   221 . log close
            name:  <unnamed>
             log:  /Users/robertyaffee/Documents/data/research/chwk/phase3/geodesicDi
      > stanceComputations.smcl
        log type:  smcl
       closed on:  20 Jan 2012, 18:05:09
      -------------------------------------------------------------------------------