function ComputeFormCD()
{
  var latdeg1,latmin1,latsgn1,londeg1,lonmin1,lonsgn1;
  var latdeg2,latmin2,latsgn2,londeg2,lonmin2,lonsgn2;
  var lat1,lat2,lon1,lon2;
  var dconv,dist,crs12,crs21;
  var argacos;
  var a,invf;
  
  /* Get coordinates of point 1 */
  latdeg1=parseFloat(document.InputFormCD.point1.value.substring(0,2));
  latmin1=parseFloat(document.InputFormCD.point1.value.substring(2,7))/1000.0;
  latsgn1=document.InputFormCD.point1.value.substring(7,8);
  lat1=(Math.PI/180)*(latdeg1+latmin1/60.0);
  if (latsgn1=='S') lat1 = -lat1;

  londeg1=parseFloat(document.InputFormCD.point1.value.substring(9,12));
  lonmin1=parseFloat(document.InputFormCD.point1.value.substring(12,17))/1000.0;
  lonsgn1=document.InputFormCD.point1.value.substring(17,18);
  lon1=(Math.PI/180)*(londeg1+lonmin1/60.0);
  if (lonsgn1=='W') lon1 = -lon1;

  /* Get coordinates of point 2 */
  latdeg2=parseFloat(document.InputFormCD.point2.value.substring(0,2));
  latmin2=parseFloat(document.InputFormCD.point2.value.substring(2,7))/1000.0;
  latsgn2=document.InputFormCD.point2.value.substring(7,8);
  lat2=(Math.PI/180)*(latdeg2+latmin2/60.0);
  if (latsgn2=='S') lat2 = -lat2;

  londeg2=parseFloat(document.InputFormCD.point2.value.substring(9,12));
  lonmin2=parseFloat(document.InputFormCD.point2.value.substring(12,17))/1000.0;
  lonsgn2=document.InputFormCD.point2.value.substring(17,18);
  lon2=(Math.PI/180)*(londeg2+lonmin2/60.0);
  if (lonsgn2=='W') lon2 = -lon2;

  /* Calculate distance and out/return courses */
  dist  = wgs84_distance(lat1,lon1,lat2,lon2);
  crs12 = ortho_bear(lat1,lon1,lat2,lon2);
  crs21 = ortho_bear(lat2,lon2,lat1,lon1);

  /* Write output to form */
  document.OutputFormCD.dist.value =format(dist,1);
  document.OutputFormCD.crs12.value=Math.round(crs12);
  document.OutputFormCD.crs21.value=Math.round(crs21);
}

function wgs84_distance(phi1, lambda1, phi2, lambda2)
{
  /* phi    = geodetic latitude  */
  /* lambda = geodetic longitude */
  
  var a = 6378.137;             /* WGS84 equatorial radius    */
  var f = 1.0/298.257223563;    /* WGS84 ellipsoid flattening */
  var EPSILON = 0.00000000005;
  var MAXITER = 100;
   
  var iter=0;
  var lambda1, phi1, lambda2, phi2 ;
  var tU1, tU2, cU1, cU2, sU1, sU2;
  var cU12, sU12, sc12, cs12;
  var L0, L, C, deltaL;
  var sL, cL, sinS, cosS, tanS, Sigma, sinA, cos2A, cosM;
  var u2, A, B, dSigma, distance;
  
  /* Reduced latitudes U1 and U2 */
  tU1 = (1.0-f)*Math.tan(phi1);         /* tan(U1) */
  tU2 = (1.0-f)*Math.tan(phi2);
  cU1 = 1.0/Math.sqrt(1.0 + tU1*tU1);   /* cos(U1) */
  cU2 = 1.0/Math.sqrt(1.0 + tU2*tU2);
  sU1 = tU1*cU1;                        /* sin(U1) */
  sU2 = tU2*cU2;
  
  /* Products to be calculated only once */
  cU12 = cU1*cU2;
  sU12 = sU1*sU2;
  cs12 = cU1*sU2;
  sc12 = sU1*cU2;
  
  /* Initial estimate for longitude difference on sphere */
  L0 = lambda2 - lambda1;
  
  /* Loop until no significant change in L */
  
  L  = L0 + 1.0;  /* force one pass */  
  while ( (Math.abs(L-L0) > EPSILON) && (iter < MAXITER) )
  {
    L0 = L;
    iter++;
    
    sL = Math.sin( L );
    cL = Math.cos( L );
    
    sinS  = Math.sqrt( Math.pow(cU2*sL,2) + Math.pow(cs12-sc12*cL,2) );
    cosS  = sU12 + cU12*cL;
    tanS  = sinS/cosS;
    Sigma = Math.acos(cosS);
    
    sinA  = cU12 * sL / sinS;
    cos2A = 1-sinA*sinA;
    cosM = cosS - (2*sU12/cos2A);
    
    C = (0.625*f)*cos2A*(4+f*(4-3*cos2A));
    L = (lambda2-lambda1) 
      + (1-C)*f*sinA*(Sigma+C*sinS*(cosM+C*cosS*(-1+2*cosM*cosM)));
  }
  
  u2 = cos2A * (1/Math.pow(1-f,2) - 1);
  A  = 1 + u2/256 * (64+u2*(-12+5*u2));
  B  = u2/512 * (128+u2*(-64+37*u2));
  dSigma = B*sinS*(cosM+0.25*B*cosS*(-1+2*cosM*cosM));
  
  distance = (1-f)*a * A * (Sigma-dSigma);
  return distance;
}


function ortho_bear(phi1, lambda1, phi2, lambda2) 
{
  /* phi    = geodetic latitude   */
  /* lambda = geodetic longitude  */

  var cosD, D ;
  var cosB, bearing ;
  
  cosD = Math.cos(phi1)*Math.cos(phi2)*Math.cos(lambda2-lambda1) + Math.sin(phi1)*Math.sin(phi2) ;
  if (cosD > 1.0) cosD = 1.0 ;
  D = Math.acos(cosD) ;
  
  if (D == 0.0) return 0.0;  /* point 1 and 2 coincide */
  
  cosB = (Math.sin(phi2)-Math.sin(phi1)*cosD)/(Math.sin(D)*Math.cos(phi1));
  if (cosB >  1.0) cosB =  1.0 ;
  if (cosB < -1.0) cosB = -1.0 ;
    
  if ( Math.sin(lambda2-lambda1) > 0.0 )
    bearing = Math.acos(cosB) ;
  else
    bearing = 2.0*Math.PI - Math.acos(cosB) ;
   
  bearing *= 180.0/Math.PI ;
  return bearing ;
}

function format(expr, decplaces)
{
  var str = ""+Math.round(eval(expr)*Math.pow(10,decplaces))
  while (str.length <= decplaces) {
    str= "0" + str
  }
  var decpoint = str.length-decplaces
  return str.substring(0,decpoint) + "." + str.substring(decpoint,str.length)
}