Great Circle Information

[toml_12953]

Another program by Stephen R. Schmitt!

Code:

distance();
course();
latitude();
longitude();

export GreatCircle()

begin

    local sgn, deg, minn, secc;
    local lat1, lat2, lat3, lat4;
    local lon1, lon2, lon3, lon4;
    local dist, crse;

    // Sydney, AS     33°53' S  151°13' E
    lat1 := latitude(  -1,  33, 53, 0 );
    lon1 := longitude( 1, 151, 13, 0 );
    
    // London, UK,    51°30' N    0°07' W
    lat2 := latitude(  1, 51, 30, 0 );
    lon2 := longitude( -1,  0,  7, 0 );

    // Boston, MA US  42°21' N   71°04' W
    lat3 := latitude(  1, 42, 21, 0 );
    lon3 := longitude( -1, 71,  4, 0 );

    // Tokyo, JA,     35°41' N  139°45' E
    lat4 := latitude(  1,  35, 41, 0 );
    lon4 := longitude( 1, 139, 45, 0 );

    print();
    print("Sydney, AS to London, UK");
    dist := distance(lat1, lon1, lat2, lon2);
    crse := course(lat1, lon1, lat2, lon2);
    print("distance = "+ round(dist,2)+ "  kilometers");
    print("heading  = "+ round(crse,2)+ "° true (initially)\n");

    print("London, UK to Boston, MA US");
    dist := distance(lat2, lon2, lat3, lon3);
    crse := course(lat2, lon2, lat3, lon3);
    print("distance = "+ round(dist,2)+ "  kilometers");
    print("heading  = "+ round(crse,2)+ "° true (initially)\n");
    
    print("Boston, MA US to Tokyo, JA");
    dist := distance(lat3, lon3, lat4, lon4);
    crse := course(lat3, lon3, lat4, lon4);
    print("distance = "+ round(dist,2)+ "  kilometers");
    print("heading  = "+ round(crse,2)+ "° true (initially)\n");

    print("Tokyo, JA to Sydney, AS");
    dist := distance(lat4, lon4, lat1, lon1);
    crse := course(lat4, lon4, lat1, lon1);
    print("distance = "+ round(dist,2)+ "  kilometers");
    print("heading  = "+ round(crse,2)+ "° true (initially)");

end;

// distance using Meeus approximation
distance( lat1, lon1, lat2, lon2 )

begin

    local F := (lat1 + lat2)/2.0;
    local G := (lat1 - lat2)/2.0;
    local L := (lon1 - lon2)/2.0;

    local sinG2 := sin(G)^2;
    local cosG2 := cos(G)^2;
    local sinF2 := sin(F)^2;
    local cosF2 := cos(F)^2;
    local sinL2 := sin(L)^2;
    local cosL2 := cos(L)^2;

    local S := sinG2*cosL2 + cosF2*sinL2;
    local C := cosG2*cosL2 + sinF2*sinL2;

    local w := atan(sqrt(S/C));
    local R := sqrt(S*C)/w;

    local a := 6378.137;            // WGS-84 equatorial radius
    local f := 1.0/298.257223563;   // WGS-84 ellipsoid flattening factor

    local D := 2*w*a;
    local H1 := (3*R - 1)/(2*C);
    local H2 := (3*R + 1)/(2*S);

    local dist := D*(1 + f*H1*sinF2*cosG2 - f*H2*cosF2*sinG2);

    return round(100*dist,0)/100.0;

end;

// course using spherical trigonometry
// does not work if one latitude is polar!!!
course( lat1, lon1, lat2, lon2 )

begin
    local L, D, cosD, C, cosC;

    L    := lon2 - lon1;
    cosD := sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(L);  
    D    := acos(cosD);
    cosC := (sin(lat2) - cosD*sin(lat1))/(sin(D)*cos(lat1));
    
    // numerical error can result in |cosC| slightly > 1.0 
    if cosC > 1.0  then
        cosC := 1.0;
    end;
    if cosC < -1.0 then
        cosC := -1.0;
    end;
 
    C := 180.0*acos(cosC)/PI;

    if sin(L) < 0.0 then
        C := 360.0 - C;
    end;

    return round(100*C,0)/100.0;

end;

// convert to radians
latitude( sgn, deg, minn, secc )

begin
    local lat := 0;

    if (0.0 <= secc) and (secc < 60.0) then
        lat := lat + secc/3600.0;
    else
        print("check latitude seconds!");
    end;

    if (0.0 <= minn) and (minn < 60.0) then
        lat := lat + minn/60.0;
    else
        print("check latitude minutes!");
    end;

    if (0.0 <= deg) and (deg <= 90.0) then
        lat := lat + deg;
    else
        print("check latitude degrees!");
    end;

    if (0.0 <= lat) and (lat <= 90.0) then
        lat := PI*lat/180.0;
    else
        print("latitude range error!");
    end;

    // assert (sgn = +1) or (sgn = -1)
    return sgn*lat;

end;

// convert to radians
longitude( sgn, deg, minn, secc )

begin
    local lon := 0;

    if (0.0 <= secc) and (secc < 60.0) then
        lon := lon + secc / 3600.0;
    else
        print("check longitude seconds!");
    end;

    if (0.0 <= minn) and (minn < 60.0) then
        lon := lon + minn / 60.0;
    else
        print("check longitude minutes!");
    end;

    if (0.0 <= deg) and (deg <= 180.0) then
        lon := lon + deg;
    else
        print("check longitude degrees!");
    end;

    if (0.0 <= lon) and (lon <= 180.0) then
        lon := PI * lon / 180.0;
    else
        print("longitude range error!");
    end;

    // assert (sgn = +1) or (sgn = -1)
    return sgn*lon;

end;