Moon Distance

#include <stdio.h>
#include <math.h>

const double pi = 2*acos(0);
const double R = 6371; /* the radius of earth in km */

/* given latitudes on same longitudes and zenith angles, get the distance of the moon.
   only works when latitudes are on different hemispheres.
   latitude on southern hemisphere is negative */
double get_moon_distance(double lat1, double lat2, double zen1, double zen2);

int main(){
    double lat1, lat2, zen1, zen2;

    while (scanf("%lf%lf%lf%lf", &lat1, &lat2, &zen1, &zen2) != EOF ){
        printf("%lf\n", get_moon_distance(lat1, lat2, zen1, zen2));
    }

    return 0;
}

double get_moon_distance(double lat1, double lat2, double zen1, double zen2) {
    double tmp;
    double lin_dist;
    double ang_dist;
    double lunar_angle_between_points;
    double moon_dist_from_center;
    double moon_dist_from_upper_point;

    /* make sure lat1 > lat2 */
    if ( lat1 < lat2 ) {
        tmp = lat1; lat1 = lat2; lat2 = tmp;
        tmp = zen1; zen1 = zen2; zen2 = tmp;
    }

    /* convert angles into radian */
    lat1 *= pi / 180.; lat2 *= pi / 180.; zen1 *= pi / 180.; zen2 *= pi / 180.;

    ang_dist = lat1 - lat2;
    lin_dist = R * sin(ang_dist) / cos(ang_dist/2);
    lunar_angle_between_points = zen1 + zen2 - ang_dist;
    moon_dist_from_upper_point = lin_dist * cos(- zen2 + ang_dist/2) / sin(lunar_angle_between_points);
    tmp = R / moon_dist_from_upper_point;
    moon_dist_from_center = moon_dist_from_upper_point * sqrt(1 + tmp * tmp + 2 * tmp * cos(zen1));
    return moon_dist_from_center - R;
}