# Started on 11.09.2000
# KS Gleditsch kgleditsch@ucsd.edu
#
# This creates a list of distances between captials
# for all dyads of states with coordinates in latlong file
#
# This version: 4 June 2003

open(LATLONG,"latlong.csv") || die "Cannot open latlong.csv";
@latlong = <LATLONG>;
open(OUT,">capdist.csv") || die "Cannot write capdist.csv";
print OUT "numa,ida,numb,idb,kmdist,midist\n";


# Modfied Fortran code received by MDW
#      function arcdis (lat1,long1,lat2,long2)
#      function to compute distance between 2 points given their
#      latitudes and longitudes, and assuming that the earth is a sphere
#      real lat1,long1,lat2, long2, (pi=3.14159 is predifined)

#     compute arc differences in radians
#      ab=((90.-lat1)*pi)/180.
#      ac=((90.-lat2)*pi)/180.
#      bc=abs(((long1-long2)*pi)/180.)
#      compute arc distance
#      arcdis=acos(cos(ab)*cos(ac)+sin(ab)*sin(ac)*cos(bc))*3969.665 

#row.names CtryNr CtryNm latg latm longg longm
for($i=0;$i<@latlong;$i++){
  $stuff = $latlong[$i];
  chop($stuff);
  @info = split(/,/,$stuff);
  $numa = $info[0];
  $ida = $info[1];
  $latg = $info[2];
  $latm = $info[3];
  $longg = $info[4];
  $longm = $info[5];
  # digitize latitude
  if($latg<0){$lat_a = ($latg-($latm/60));
  } else { $lat_a = ($latg+($latm/60));}
  # digitize longitude
  if($longg<0){$long_a = ($longg-($longm/60));
  } else { $long_a = ($longg+($longm/60));}
 
  for($j=0;$j<@latlong;$j++){
    unless($i==$j){
       $stuff = $latlong[$j];
       chop($stuff);
       @info = split(/,/,$stuff);
       $numb = $info[0];
       $idb = $info[1];
       $latg = $info[2];
       $latm = $info[3];
       $longg = $info[4];
       $longm = $info[5];
       # digitize latitude
       if($latg<0){$lat_b = ($latg-($latm/60));
       } else { $lat_b = ($latg+($latm/60));}
       # digitize longitude
       if($longg<0){$long_b = ($longg-($longm/60));
       } else { $long_b = ($longg+($longm/60));}


       $pi = atan2(1,1) * 4;
       $km_constant = 6367;
       $miles_constant = 3970;

       sub acos{ atan2( sqrt(1-$_[0] * $_[0]), $_[0] ) }

       $ab = ((90-$lat_a)*$pi)/180;
       $ac = ((90-$lat_b)*$pi)/180;
       $bc = abs((($long_a-$long_b)*$pi)/180);
       $km_distance =  sprintf "%.0f", (acos(cos($ab)*cos($ac)+sin($ab)*sin($ac)*cos($bc))*$km_constant);
       $miles_distance =  sprintf "%.0f", (acos(cos($ab)*cos($ac)+sin($ab)*sin($ac)*cos($bc))*$miles_constant);

       print OUT join(",",($numa,$ida,$numb,$idb,$km_distance,$miles_distance))."\n";


    } # End unless
   } #end loop over states


}