*
* The first order differential equation of the smooth geodetic latitude
* as a function of longitude is solved by moving from one point to the
* next on a finite grid of equidistant sampling points on the longitude,
* using a predictor-only Newton approximation.
* The parameter which defines the initial direction in the local
* tangential plane is gathered by a shooting method (a starting
* guess is obtained from a spherical approximation) and a fixed number of
* iteration through the parameter space.
*
* @since 2008-04-12
* @see preprint Geodetic line at constant altitude above the Ellipsoid,
* arXiv:0711.0642 [math.MG]
* @author Richard J. Mathar
*/
class Geod
{
/** The inverse flattening factor of the IERS 2003 convention.
*/
static public final double IERS_TN32_FLAT= 298.25642 ;
/** The inverse flattening factor of the IERS 1996 convention.
*/
static public final double IERS_TN21_FLAT= 298.25642 ;
/** The inverse flattening factor of the GRS80 model.
*/
static public final double GRS80_FLAT= 298.257222101 ;
/** The inverse flattening factor of the WGS84 model.
*/
static final double WGS84_FLAT= 298.257223563 ;
/** Equatorial radius of the IERS 2003 convention, meters.
*/
static public final double IERS_TN32_RHO_E=6378136.6 ;
/** Equatorial radius of the IERS 1996 convention, meters.
*/
static public final double IERS_TN21_RHO_E=6378136.49 ;
/** Equatorial radius of the GRS 1980 system, meters.
*/
static public final double GRS80_RHO_E=6378137.0 ;
/** Equatorial radius of the WGS84, meters.
*/
static public final double WGS84_RHO_E=6378137.0 ;
/** Three enumeration values to convert angular units on user input and output
*/
interface AngleUnit {
int RAD =0, DEG = 1, GON = 2;
} ;
/** Equatorial radius in meters.
*/
double rho_equat ;
/** Eccentricity.
* For the Earth close to 0.08.
*/
double eccen ;
/** Altitude of the inverse problem above the ellipsoid.
* Value equals 0 on the ellipsoid surface. May have both signs.
* Measured in the same units as Geod::rho_equat.
*/
double altit ;
/** Default Constructor.
* The variables are set to the defaults of the surface of the WGS84.
*/
public Geod()
{
this(WGS84_RHO_E,Math.sqrt((2.0-1./WGS84_FLAT)/WGS84_FLAT),0.) ;
}
/** Constructor.
* @param rho equatorial radius [m]
* @param e eccentricity
* @param h altitude above the reference ellipsoid [m]
* @todo catch error conditions of e>=1, rho<0 or h< polar radius.
*/
public Geod(double rho, double e, double h)
{
rho_equat = rho ;
eccen = e ;
altit = h ;
}
/** Convert geodetic to Cartesian coordinates.
* @param tau the sine of the geodetic latitude
* @param lambd the geodetic longitude [rad]
* @return the three components of the Cartesian coordinates
* in the same units as rho_equat and altit.
* @see eq (15) of the preprint
*/
public double[] getCartesian(final double tau, final double lambd)
{
double[] cart = new double[3] ;
final double N = curvN(tau) ;
final double sinphi = Math.sqrt(1.-tau*tau) ;
cart[0] = (N+altit)*sinphi*Math.cos(lambd) ;
cart[1] = (N+altit)*sinphi*Math.sin(lambd) ;
cart[2] = (N*(1.+eccen)*(1.-eccen)+altit)*tau ;
return cart ;
}
/** Curvature parameter N [m]
* @param tau sine of the latitude
* @return N
according to eq (13) of the preprint.
*/
public double curvN(final double tau)
{
return rho_equat/Math.sqrt( (1.+eccen*tau)*(1.-eccen*tau) ) ;
}
/** Derivative d N/d tau [m].
* @param tau sine of the latitude
* @return the derivative of N(tau)
*/
double dNdtau(double tau)
{
return curvN(tau)*eccen*eccen*tau/((1.+eccen*tau)*(1.-eccen*tau)) ;
}
/** Convert eccentricity to flattening factor.
* @param e the (first) eccentricity
* @return f = 1-sqrt(1-eccentricity^2).
*/
public static double flatt(final double e)
{
final double e2 = e*e ;
if ( e2 < 3.8e-3)
/* a Taylor approximation up to O(e^14), relative accuracy to E-16. */
return e2*(0.5+e2*(1./8.+e2*(1./16.*e2*(5./128.+e2*(7./256.+e2*(21./1024.+e2*33./2048.)))))) ;
else
return 1.-Math.sqrt(1.-e2) ;
}
/** Flattening factor.
* @return f = 1-sqrt(1-eccentricity^2).
*/
public double flatt()
{
return flatt(eccen) ;
}
/** Equatorial radius of curvature M [m]
* @param tau sine of the latitude
* @return M
according to eq (18) of the preprint.
*/
public double curvM(final double tau)
{
return curvN(tau)*(1.-eccen*eccen)/( (1.+eccen*tau)*(1.-eccen*tau) ) ;
}
/** Derivative d M /d tau [m]
* @param tau sine of the latitude
* @return the derivative of curvM() with respect to the parameter tau.
*/
public double dMdtau(final double tau)
{
return 3.*curvM(tau)*eccen*eccen*tau/( (1.+eccen*tau)*(1.-eccen*tau) ) ;
}
/** second Derivative d^2 M /d tau^2 [m]
* @param tau sine of the latitude
* @return the second derivative of curvM() with respect to the parameter tau.
*/
public double d2Mdtau2(final double tau)
{
final double etau2 = eccen*eccen*tau*tau ;
return 3.*curvM(tau)*Math.pow(eccen/(1.-etau2),2.)*(1.+4.*etau2) ;
}
/** Gauss parameter E [m^2].
* @param tau sine of the latitude
* @return E
according to eq (19) of the preprint.
*/
double GaussE(final double tau)
{
return Math.pow(curvN(tau)+altit,2.)*(1.-tau*tau) ;
}
/** Derivative d E/ d tau [m^2]
* @param tau sine of the latitude
* @return the derivative of GaussE() with respect to the parameter tau
*/
public double dEdtau(final double tau)
{
return -2.*tau*(curvN(tau)+altit)*(curvM(tau)+altit) ;
}
/** Derivative d^2 E/ d tau^2 [m^2]
* @param tau sine of the latitude
* @return the 2nd derivative of GaussE() with respect to the parameter tau
*/
public double d2Edtau2(final double tau)
{
/* Use the Leibniz rule applied to the first derivative, which was -2*tau*(N+h)(M+h)
*/
final double Nh = curvN(tau)+altit ;
final double Mh = curvM(tau)+altit ;
return -2.*(Nh*Mh+tau*dNdtau(tau)*Mh+tau*Nh*dMdtau(tau)) ;
}
/** Gauss parameter G [m^2].
* @param tau sine of the latitude
* @return G
according to eq (21) of the preprint.
* @since 2008-04-22
*/
double GaussG(final double tau)
{
return Math.pow(curvM(tau)+altit,2.) ;
}
/** Minimum polar distance.
* @param c3 the constant c3 along one individual trajectory
* @return the value of tau_m
according to eq (66) of the preprint
that could be reached if the trajectory were
* closed over the surface. That is the maximum of |sin(phi)|, or the value
* of |sin(phi)| where the derivative of sin(phi) with respect to the longitude
* becomes zero.
* @since 2008-04-22
*/
double tauTropic(final double c3)
{
final double c3sqr = c3*c3 ;
final double H = rho_equat+altit ;
final double rhoebar = rho_equat/H ;
final double esqr = eccen*eccen ;
/* Initial estimate from Taylor expansion in orders of squared eccentricity
* See eq (67) of the preprint
*/
double tau2 = (1.-c3sqr/(H*H))*(1+c3sqr*rhoebar*esqr/(H*H)*
(1.+0.25*(3.*(1.-c3sqr/(H*H))+rhoebar*(7.*c3sqr/(H*H)-3.))*esqr)
) ;
final double hsqr = altit*altit ;
final double rhosqr = rho_equat*rho_equat ;
/* coefficients in the 4th order polynomial of tau^2 which is to be solved */
double[] coeff = new double[5] ;
coeff[0] = (Math.pow(rho_equat+altit,2.)-c3sqr)*(Math.pow(rho_equat-altit,2.)-c3sqr) ;
coeff[1] = 2.*(-Math.pow(rhosqr-hsqr,2.)
+c3sqr*(rhosqr+hsqr)
+esqr*(-Math.pow(c3sqr-hsqr,2.)+rhosqr*(c3sqr+hsqr))
) ;
coeff[2] = Math.pow(rhosqr-hsqr,2.)
+2.*esqr*(2.*hsqr*hsqr-2.*c3sqr*hsqr-2.*rhosqr*hsqr-rhosqr*c3sqr)
+esqr*esqr*Math.pow(c3sqr-hsqr,2.) ;
coeff[3] = 2.* Math.pow(altit*eccen,2.)*(rhosqr-hsqr+esqr*(c3sqr-hsqr)) ;
coeff[4] = Math.pow(altit*eccen,4.) ;
/* Small number of loops for self-consistency, using some Horner scheme for
* the function to be zeroed and its first derivative */
for(int loop=0; loop < 4 ;loop++)
{
System.out.println("# taum2 "+tau2) ;
tau2 -= ( coeff[0]+tau2*(coeff[1]+tau2*(coeff[2]+tau2*(coeff[3]+tau2*coeff[4])))
)/(
coeff[1]+tau2*(2.*coeff[2]+tau2*(3.*coeff[3]+tau2*4.*coeff[4]))
) ;
}
return Math.sqrt(tau2) ;
}
/** Derivative d tau/d lambda.
* @param tau sine of the latitude
* @param c3 the constant c3 along one individual trajectory
* @return absolute value of the derivative
according to eq (64) of the preprint.
*/
double dtaudlambda(double tau, double c3)
{
/* treat the factors (N+h)^2*(1-tau^2) and the sqrt(..)/(h+M) separately */
return GaussE(tau)*dtauds(tau,c3)/c3 ;
}
/** Derivative d s/d lambda.
* @param tau sine of the latitude
* @param c3 the constant c3 along one individual trajectory
* @return absolute value of the derivative of the path length
according to eq (63) of the preprint.
*/
double dsdlambda(double tau, double c3)
{
return GaussE(tau)/c3 ;
}
/** Discriminant under the square root of d tau / ds.
* @param tau sine of the latitude
* @param c3 the constant c3 along one individual trajectory
* @return 1-tau^2-(c3/(N+h))^2.
*/
double discrT(final double tau, final double c3)
{
/* N+altitude */
final double Nh = curvN(tau)+altit ;
return 1.-tau*tau-Math.pow(c3/Nh,2.) ;
}
/** Derivative of T with respect to tau.
* @param tau sine of the latitude
* @param c3 the constant c3 along one individual trajectory
* @return derivative of discrT() with respect to tau
*/
double dTdtau(final double tau, final double c3)
{
final double N = curvN(tau) ;
final double Nh = N+altit ;
return 2*tau*( N*Math.pow(c3*eccen/Nh,2.)/(Nh*(1.+eccen*tau)*(1.-eccen*tau)) -1. ) ;
}
/** Second derivative of T with respect to tau.
* @param tau sine of the latitude
* @param c3 the constant c3 along one individual trajectory
* @return the second derivative of discrT() with respect to tau
*/
double d2Tdtau2(final double tau, final double c3)
{
final double N = curvN(tau) ;
final double Nh = N+altit ;
final double e2 = eccen*eccen ;
final double Oneetau = (1.+eccen*tau)*(1.-eccen*tau) ;
return -2.*( 1.-N*c3*c3*e2*(1.+3.*altit*e2*tau*tau/Nh/Oneetau)/Math.pow(Nh,3.)/Oneetau ) ;
}
/** Derivative d tau/d s [1/m].
* @param tau sine of the latitude
* @param c3 the constant c3 along one individual trajectory
* @return absolute value of the derivative
according to eq (61) of the preprint.
*/
double dtauds(double tau, double c3)
{
return Math.sqrt( discrT(tau,c3) )/(altit+curvM(tau)) ;
}
/** Second derivative d^2 tau/d lambda^2.
* @param tau sine of the latitude
* @param c3 the constant c3 along one individual trajectory
* @return value of the derivative
of eq (64) of the preprint.
*/
double d2taudlambda2(double tau, double c3)
{
/* The array contains values of E, sqrt(T) and 1/(h+M) in the components facto[][0],
* and their first derivatives with respect to tau in facto[][1].
*/
double[][] facto = new double[3][2] ;
facto[0][0] = GaussE(tau) ;
facto[1][0] = Math.sqrt(discrT(tau,c3) ) ;
facto[2][0] = 1./( curvM(tau)+altit) ;
facto[0][1] = dEdtau(tau) ;
/* derivative of square root is one half of the inverse multiplied by interior derivative
*/
facto[1][1] = dTdtau(tau,c3)/(2.*facto[1][0]) ;
/* derivative of 1/(h+M) is negative of its square multiplied by interior derivative
*/
facto[2][1] = -dMdtau(tau)*facto[2][0]*facto[2][0] ;
/* The Leibniz rule of derivatives is applied to collect the derivative of
* the product of facto[i][0], i=0..2.
*/
final double resul = facto[0][1]*facto[1][0]*facto[2][0] +facto[0][0]*facto[1][1]*facto[2][0]
+facto[0][0]*facto[1][0]*facto[2][1] ;
/* resul contains so far the derivative (d/dtau) (d tau/d lambda), without
* the c3 in the denominator. Postmultiply with d tau/ d lambda to generate
* the d^2 tau/ d lambda^2.
*/
return dtaudlambda(tau,c3)*resul/c3 ;
}
/** Third derivative d^3 tau/d lambda^3.
* @param tau sine of the latitude
* @param c3 the constant c3 along one individual trajectory
* @return value of the 2nd derivative
of eq (64) of the preprint.
*/
double d3taudlambda3(double tau, double c3)
{
final double T = discrT(tau,c3) ;
/* Faa di Bruno formula:
* d^3 tau/d lambda^3 = (d^2tau/d lambda^2)* (d/dtau) (dtau/dlambda)+ (dtau/dlambda)^2* (d^2/dtau^2) (dtau/dlambda
*/
/* The array contains values of E, sqrt(T) and 1/(h+M) in the components [][0],
* and their first derivatives with respect to tau in [][1], 2nd derivatives in [][2].
*/
double[][] facto = new double[3][3] ;
facto[0][0] = GaussE(tau) ;
facto[1][0] = Math.sqrt(T) ;
facto[2][0] = 1./( curvM(tau)+altit) ;
facto[0][1] = dEdtau(tau) ;
/* derivative of square root is one half of the inverse multiplied by interior derivative
*/
final double Tprime = dTdtau(tau,c3) ;
facto[1][1] = Tprime/(2.*facto[1][0]) ;
/* derivative of 1/(h+M) is negative of its square multiplied by interior derivative
*/
final double Mprime = dMdtau(tau) ;
facto[2][1] = -Mprime*facto[2][0]*facto[2][0] ;
facto[0][2] = d2Edtau2(tau) ;
/* First derivative of sqrt(T) was dTdtau/(2*sqrt(T)). The 2nd derivative therefore is
* [ (d^2/dtau^2) T -(d T/d tau)^2/(2T)]/[2sqrt(T)], which we store in facto[1][2].
*/
facto[1][2] = ( d2Tdtau2(tau,c3) -Tprime*Tprime/( 2.*T ) )/(2.*facto[1][0]) ;
/* Second derivative of 1/(h+M). First one was -(dM/dtau)/(h+M)^2. Second one is
/* -(d^2M/dtau^2)/(h+M)^2+2(dM/dtau)^2/(h+M)^3.
*/
facto[2][2] = ( 2.*Mprime*Mprime*facto[2][0]-d2Mdtau2(tau) )*facto[2][0]*facto[2][0] ;
/* The value of (d/d tau) (d tau/d lambda). The implementation is not
* tuned for quick evaluation, since the same value is calculated again in #d2taudlambda2.
*/
final double dtaudlambdadtau = ( facto[0][1]*facto[1][0]*facto[2][0] +facto[0][0]*facto[1][1]*facto[2][0]
+facto[0][0]*facto[1][0]*facto[2][1] )/c3 ;
/* the value of (d^2/dtau^2) (dtau/dlambda) with Leibniz rule, multinomial case
*/
final double d3taudlambdadtau2 = ( facto[0][2]*facto[1][0]*facto[2][0] +2.*facto[0][1]*facto[1][1]*facto[2][0]
+2.*facto[0][1]*facto[1][0]*facto[2][1] +facto[0][0]*facto[1][2]*facto[2][0]
+2.*facto[0][0]*facto[1][1]*facto[2][1] +facto[0][0]*facto[1][0]*facto[2][2] )/c3 ;
/* Faa di Bruno rule of derivatives is applied to collect the derivative
* (d^2/dlambda^2) dtau/dlambda
* = (d^2 tau/dlambda^2)*(d/dtau) (d tau/d lambda)+(d tau/d lambda)^2* (d^2/d tau^2) (d tau/d lambda)
*/
return d2taudlambda2(tau,c3)*dtaudlambdadtau + Math.pow(dtaudlambda(tau,c3),2.)* d3taudlambdadtau2;
}
/** Second derivative d^2 s/d lambda^2 = (d/d lambda) (E/c3).
* @param tau sine of the latitude
* @param c3 the constant c3 along one individual trajectory
* @return value of the 2nd derivative
of eq (63) of the preprint.
*/
double d2sdlambda2(double tau, double c3)
{
/* First the derivative (d/dtau) (d s/d lambda), without
* the c3 in the denominator. Postmultiply with d tau/ d lambda to generate
* the d^2 s/ d lambda^2.
*/
return dtaudlambda(tau,c3)*dEdtau(tau)/c3 ;
}
/** Third derivative d^3 s /d lambda ^3.
* @param tau sine of the latitude
* @param c3 the constant c3 along one individual trajectory
* @return value of the third derivative
* @since 2010-11-26
of eq (67) of the preprint
*/
double d3sdlambda3(double tau, double c3)
{
final double dt2 = dtaudlambda(tau,c3) ;
return ( d2taudlambda2(tau,c3)*dEdtau(tau) + Math.pow(dt2,2)*d2Edtau2(tau) )/c3 ;
}
/** Spherical approximation to parameter c3
* @param phi the latitudes at start and end [rad]
* @param lambd the longitudes at start and end [rad]
* @return the parameter c3 obtained form the spherical approximation.
* @see eq (A19) in the preprint.
*/
double c3Sphere(final double[] phi, final double lambd[])
{
/* cosine of the angle Z */
final double cosZ = Math.sin(phi[0])*Math.sin(phi[1])
+Math.cos(phi[0])*Math.cos(phi[1])*Math.cos(lambd[1]-lambd[0]) ;
return (altit+rho_equat)*Math.cos(phi[0])*Math.cos(phi[1])*Math.sin(lambd[1]-lambd[0])/Math.sqrt(1.-cosZ*cosZ) ;
}
/** Approximate sign d tau / ds
* @param phi the latitudes at start and end [rad]
* @param lambd the longitudes at start and end [rad]
* @return the sign of the square root
in eq (64) of the preprint,
as obtained from the spherical approximation.
*/
double dtaudsSignum(final double[] phi, final double lambd[])
{
/* cosine of the angle Z */
final double cosZ = Math.sin(phi[0])*Math.sin(phi[1])
+Math.cos(phi[0])*Math.cos(phi[1])*Math.cos(lambd[1]-lambd[0]) ;
return Math.signum(Math.sin(phi[1])-Math.sin(phi[0])*cosZ) ;
}
/** Adjust longitude of the final point.
* @param lambd the longitudes at start and end [rad]
* The value of lambd[1] is adjusted to the range lambd[0]+-pi to ensure that the
* trajectory chooses a path along the correct hemisphere.
*/
static void adjLambdaEnd(double lambd[])
{
final double dl = lambd[1]-lambd[0] ;
if ( Math.abs( dl+2.*Math.PI) < Math.abs(dl) )
lambd[1] += 2.*Math.PI ;
if ( Math.abs( dl-2.*Math.PI) < Math.abs(dl) )
lambd[1] -= 2.*Math.PI ;
}
/** Nautical course in the local oblique horizontal [rad]
* @param tau sine of the latitude
* @param c3 the directional parameter of the solution to the inverse problem [m]
* @param north sign of the square root, +1 or -1
* @return the direction angle (North over East) in the local oblique tangential plane [rad]
* @see eq (78) of the preprint
*/
double nautAngle(final double tau, final double c3, final double north)
{
/* sinkappa and coskappa are the sine and cosine of the angle multiplied by
* a common positive factor, which retains the correct quadrant of the solution.
*/
final double sinkappa = c3/(curvN(tau)+altit) ;
final double coskappa = north* Math.sqrt( discrT(tau,c3) ) ;
return Math.atan2(sinkappa,coskappa) ;
}
/** Compute one trajectory over the lambda interval, assuming c3 given.
* @param phi start and (target) value of the geodetic latitude [rad]
* @param lambd start and end value of the longitude [rad]
* @param c3 fixed parameter of c3
* @param Nsampl number of steps into which the interval lambd[0]-lambd[1] is divided
* @param usampl if positive, each usampl'th point of the result is printed
* @param useRad one value of the AngleUnit to indicate which units are preferred for angles
* @param taylOrd the order of the Taylor approximation in the FEM.
* A value >=3 triggers use of a third order approximation, all others use of a 2nd order approximation.
* @return the overshooting of the value of tau versus the sine of phi[1] [rad]
*/
double tauShoot(final double[] phi, final double lambd[], final double c3, final int Nsampl, final int usampl,
final int useRad, final int taylOrd)
{
/* step width in the longitudes, including a sign */
final double dlam = (lambd[1]-lambd[0])/Nsampl ;
/* current value of tau = sin(phi), initialized with the value at the start */
double tau = Math.sin(phi[0]) ;
/* current value of s, the length along the line, initialized with the value at the start */
double s = 0. ;
/* the sign which chooses a N or a S route. +1 or -1 */
double north = dtaudsSignum(phi,lambd) ;
if ( usampl > 0 )
System.out.println("\n\n") ;
double kappa = nautAngle(tau,c3,north) ;
System.out.println("# start course "+ kappa +" rad, " + Math.toDegrees(kappa)+ " deg") ;
/* loop over the finite elements */
for(int i=0; i < Nsampl; i++)
{
/* in each of the elements, a 2nd or 3rd order Taylor approximation
* for dtau/ds is applied
with (64) of the preprint
*/
double deriv1 = north*dtaudlambda(tau,c3) ;
/* below, the 'north' factor appears twice and cancels */
double deriv2 = d2taudlambda2(tau,c3) ;
double deriv3 = (taylOrd >= 3) ? north*d3taudlambda3(tau,c3) : 0.;
/* Taylor extrapolation. Current value +(df/dx)*delta+(1/2)(d^2/dx^2)*delta^2+(1/6)*(d^3/dx^3)*delta^3
* Update to the value at the right end of the interval.
*/
tau += dlam*(deriv1+dlam*(0.5*deriv2+dlam*deriv3/6.) ) ;
/* We may have passed by a maximum of the equatorial distance
* and sense this by comparing the location of the extremum on the current parabolic
* estimator for tau(lambda) with the current interval. Setting the derivative of the
* Taylor estimator to zero, (df/dx)+(d^2/dx^2)*(lambda-current lambda)=0, the
* location of the extremum is -(df/dx)/( d^2f/dx^2) relative to the current lambda.
* @todo implement the 3rd order Taylor Approximation as well.
*/
final double lamFlat = -deriv1/deriv2 ;
/* in each of the elements, a 2nd or 3rd order Taylor approximation
* for ds/dlambda is applied
with (64) of the preprint
*/
double deriv1s = dsdlambda(tau,c3) ;
double deriv2s = d2sdlambda2(tau,c3) ;
double deriv3s = (taylOrd >= 3) ? d3sdlambda3(tau,c3) : 0.;
/* Update of the path length along the trajectory in the Taylor approximation
*/
s += dlam*(deriv1s+dlam*(0.5*deriv2s+dlam*deriv3s/6.) ) ;
if ( usampl > 0 )
if ( i % usampl == 0 || i == Nsampl-1 )
{
/* longitude at right end of interval [rad] */
final double longi = lambd[0]+(i+1)*dlam ;
double[] cart = getCartesian(tau, longi ) ;
/* Print the x,y,z Cartesian coordinates as the first 3 columns per line */
System.out.print(""+cart[0]+" "+cart[1]+" "+cart[2]) ;
kappa = nautAngle(tau,c3,north) ;
/* Print longitude, geodetic latitude, nautical angle as columns 4 to 6 */
if ( useRad == AngleUnit.RAD )
System.out.print(" "+ longi +" "+ Math.asin(tau)+ " " +kappa ) ;
else if ( useRad == AngleUnit.DEG )
System.out.print(" "+ Math.toDegrees(longi) + " " + Math.toDegrees(Math.asin(tau))
+" "+Math.toDegrees(kappa) ) ;
else
System.out.print(" "+ 10.*Math.toDegrees(longi)/9. + " "
+ 10.*Math.toDegrees(Math.asin(tau))/9.
+" "+10.*Math.toDegrees(kappa)/9. ) ;
/* Print length of trajectory in column 7 */
System.out.print(" "+s) ;
System.out.println("") ;
}
if ( lamFlat/dlam >= 0. && lamFlat/dlam < 1.)
{
north *= -1. ;
}
}
/* tau now is the estimate of the final position; return the error */
return tau-Math.sin(phi[1]) ;
}
/** Perform 4 iterations on adjusting the parameter c3.
* @param phi the latitudes at start and end [rad]
* @param lambd the longitudes at start and end [rad]
* @param Nsampl the number of divisions of the longitude axis
* @param usampl undersampling factor for the printout in the converged trajectory
* @param useRad the unit of the angles which are reported
* @param taylOrd the order of the Taylor expansion of the core integral, 2 or 3.
* @return a convergent to the parameter c3
*/
public double c3shoot(final double[] phi, final double lambd[],final int Nsampl, final int usampl, final int useRad,
final int taylOrd)
{
/* For a parabolic (2nd order) approximation of the shooting
* error (units of tau) use 3 abscissa and target missing parameters
*/
double[] err = new double[4] ;
double[] c3spread = new double [4] ;
/* start with the estimate from the spherical approximation, zero eccentricity */
c3spread[0] = c3Sphere(phi,lambd) ;
err[0] = tauShoot(phi,lambd,c3spread[0],Nsampl,-1,useRad,taylOrd) ;
System.out.println("# c3 (spherical) "+c3spread[0]+" err "+err[0]) ;
/* A stepping parameter for the variation */
final double c3Delt = -c3spread[0]*0.005 ;
c3spread[1] = c3spread[0]+c3Delt ;
err[1] = tauShoot(phi,lambd,c3spread[1],Nsampl,-1,useRad,taylOrd) ;
/* estimate of the 3rd c3 by a linear interpolation
* between the previous 2 to get zero error. The ansatz is err=err[0]+(err[1]-err[0])/c3Delt*(c-c[0]) = 0,
* where (err[1]-err[0])/c3Delt is the slope of the error curve, which is solved for c.
*/
c3spread[2] = c3spread[0]-c3Delt*err[0]/(err[1]-err[0]) ;
err[2] = tauShoot(phi,lambd,c3spread[2],Nsampl,-1,useRad,taylOrd) ;
System.out.println("# c3 (linear extrapol) "+c3spread[2]+" err "+err[2]) ;
/* Lagrange quadratic interpolation between (c3spred[0],err[0]), (c3spread[1],err[1])
* and (c3srepad[2],err[2]): err = eofCinter[0]+eofCinter[1]*c3spread+eofCinter[2]*c3spread^2
* up to a common factor which is omitted.
*/
double[] eofCinter = new double [3] ;
eofCinter[0] = err[2]*c3spread[0]*c3spread[1]*(c3spread[0]-c3spread[1])
+err[1]*c3spread[0]*c3spread[2]*(c3spread[2]-c3spread[0])
+err[0]*c3spread[1]*c3spread[2]*(c3spread[1]-c3spread[2]) ;
eofCinter[1] = err[2]*(Math.pow(c3spread[1],2.)-Math.pow(c3spread[0],2.))
+ err[1]*(Math.pow(c3spread[0],2.)-Math.pow(c3spread[2],2.))
+ err[0]*(Math.pow(c3spread[2],2.)-Math.pow(c3spread[1],2.)) ;
eofCinter[2] = err[2]*(c3spread[0]-c3spread[1])
+err[1]*(c3spread[2]-c3spread[0])
+err[0]*(c3spread[1]-c3spread[2]) ;
/* Set this quadratic to zero; normalize to monic polynomial */
eofCinter[0] /= eofCinter[2] ;
eofCinter[1] /= eofCinter[2] ;
eofCinter[2] = Math.sqrt(0.25*Math.pow(eofCinter[1],2)-eofCinter[0]) ;
/* From the two solutions take the one closer to the previous estimate */
c3spread[3] = (c3spread[2]+0.5*eofCinter[1] > 0.) ? -0.5*eofCinter[1]+eofCinter[2] : -0.5*eofCinter[1]-eofCinter[2] ;
err[3] = tauShoot(phi,lambd,c3spread[3],Nsampl,usampl,useRad,taylOrd) ;
System.out.println("# c3 (quadratic extrapol) "+c3spread[3]+" err "+err[3]) ;
return c3spread[3] ;
}
/**
* usage:
*
* java -jar Geod.jar [-h altitude/m] [-s #longitude steps] [-u sfac] [-R equat radius/m] [-e eccent] [-r]
* [-g] [-2] phi1 lam1 phi2 lam2
*
*
* Command line parameters
*