//package hehl.gl;

/**
recursive computation of meridian arcs.
author    hehl@tfh-berlin.de
version   17-May-2005
*/
public final class MB { // german for MeridianBogen, i.e. meridian arc
  private static final int MAX = 8; // largest order that makes sense
  private Ellipsoid _ell; // the ellipsoid of revolution
  private int _n;  // current order of series expansion
  private double _k2[]; // vector of coefficients in series expansion

/**
construct meridian arc object from user-defined parameters.
@param elln name of ellipsoid
@param a semimajor axis in meters
@param f flattening (a-b)/a
@param max order of evaluation (range: 1 to 8)
@throws Msg
*/
  public MB(String elln, double a, double f, int max) throws Msg {
                init(new Ellipsoid(elln,a,f),max);
        }

/**
construct meridian arc object from user-defined parameters.
@param ell an ellipsoid
@throws Msg
*/
  public MB(Ellipsoid ell)  {
                init(ell,MAX);
        }

/**
construct meridian arc object from user-defined parameters.
@param ell an ellipsoid
@param max order of evaluation (range: 1 to 8)
@throws Msg
*/
  public MB(Ellipsoid ell, int max)  {
                init(ell,max);
        }

/**
initialize meridian arc object from user-defined parameters.
@param ell the ellipsoid
@param max order of evaluation (range: 1 to 8)
*/
  private void init(Ellipsoid ell, int max) {
   _ell = ell;
   if(max < 1) max = 1; if(max > MAX) max = MAX;
   _n = max;
   _k2 = new double [order()];
   for(int n=0;n<order();n++) _k2[n] = 0.0;
   double c = 1.0;
   for(int n=1;n<=order();n++) {
     double n2=2.*n;
     c *= (n2-1.)*(n2-3.)/n2/n2 * e2();
     for(int m=0;m<n;m++) _k2[m] += c;
   }
  }

/**
construct meridian arc object from user-defined parameters with ultimate accuracy.
@param elln name of ellipsoid
@param a semimajor axis in meters
@param f flattening (a-b)/a
@throws Msg
*/
  public MB(String elln, double a, double f) throws Msg {
    this(elln,a,f,MAX);
  }

/**
construct meridian arc object on a sphere; ultimate accuracy.
@param r radius of sphere in meters
@throws Msg
*/
  public MB(double r) throws Msg {
    this("Kugel",r,0.0,MAX);
  }

/**
construct meridian arc object on Bessel ellipsoid with ultimate accuracy.
*/
  public MB() {
    init(Ellipsoid.BESSEL,MAX);
  }

/** @return the current ellipsoid */
  public Ellipsoid ell() { return(_ell); }

/** @return the name of current ellipsoid */
  public String ellname() { return(_ell.name()); }

/** @return the order of evaluation */
  public int order() { return(_n); }

/** @return the semimajor axis of ellipsoid */
  public double a()  { return(_ell.a()); }

/** @return the flattening */
  public double f()  { return(_ell.f()); }

/** @return the square of first numerical eccentricity */
  public double e2() { return(f()*(2.-f())); }

/** @return true if sphere */
  public boolean isSphere() { return(_ell.isSphere()); }

/** @return the first coefficient of series expansion */
  public double koeff1() { return(1.+_k2[0]); }

/** @return the second coefficient of series expansion */
  public double koeff2(double beta_rad) {
    double cos2 = Math.cos(beta_rad)*Math.cos(beta_rad);
    double result = _k2[0];
    double k = 1.0;
    for(int n=1;n<order();n++) {
      k *= (2.*n)/(2.*n+1.0)*cos2;
      result += _k2[n]*k;
    }
    return(result);
  }

/** @return the complex second coefficient of series expansion */
  public Complex koeff2(Complex beta) {
    Complex cos2 = Complex.cos(beta).times(Complex.cos(beta));
    Complex result = new Complex(_k2[0],0.0);
    Complex k = new Complex(1.0,0.0);
    for(int n=1;n<order();n++) {
      k = k.times((2.*n)/(2.*n+1.0)).times(cos2);
      result = result.plus(k.times(_k2[n]));
    }
    return(result);
  }

/** @return the circumference of meridian ellipse */
  public double circumference() { return(4.*smax()); }

/** @return 0.25*circumference of meridian ellipse */
  public double smax() { return(a()*koeff1()*Math.PI/2.); }

/**
computes the meridian arc from equator to point with ellipsoidal latitude phi.
@param phi the ellipsoidal latitude
@return the meridian arc in meters
*/
  public double arc(Angle phi) {
    double beta = Math.atan((1.-f())*Math.tan(phi.rad()));
    return(a()*beta*koeff1()+a()/2.*Math.sin(2.*beta)*koeff2(beta));
  }

/**
computes the meridian arc from equator to point with ellipsoidal latitude phi.
@param phideg the ellipsoidal latitude in degrees
@return the meridian arc in meters
*/
  public double arc(double phideg) {
                double s = 0.0;
                try {
                  s = arc(new Angle(phideg,Angle.DEG));
                }
                catch(Exception ex) {}
                return s;
  }

/**
computes the ellipsoidal latitude from meridian arc length.
@param s the meridian arc length in meters
@return the ellipsoidal latitude
*/
  public Angle lat(double s) throws Msg {
    final int MAXITER=10;
    double beta0 = s/a()/koeff1();
    double beta=beta0, betaold = 1.E30;
    int iter = 0;
    while(++iter < MAXITER && Math.abs(beta-betaold)>1.E-15) {
      betaold = beta;
      beta = beta0 - koeff2(beta)/2./koeff1()*Math.sin(2.*beta);
    }
    if(iter == MAXITER)
        throw new Msg("MB.lat","divergence!  iter= "+iter);
    return(new Angle(Math.atan(Math.tan(beta)/(1.-f())),Angle.RAD));
  }

        public Complex arc(Complex b) {
                Complex result = new Complex();
                try {
                        Complex test = new Complex(Complex.tan(b).times(1.-f()));
        Complex beta = Complex.atan(test);//Complex.tan(b).times(1.-f()));
        Complex koeff2complex = koeff2(beta);
        result = beta.times(a()*koeff1());
        result = result.plus(koeff2complex.times(Complex.sin(beta.times(2.))).times(a()/2.));
                }
                catch(Msg m) {
                        System.err.println(m);
                }
                catch(Exception ex) {}
                return result;
        }

/** short description of an MB object */
  public String toString() {
    return("meridian arcs on "+ellname()+", order N="+order());
  }

        public Complex phidlam2GK(Angle phi, Angle dlam)
                        throws Msg {
                return arc(Complex.complexLatitude(ell(),phi,dlam));
        }

        public Complex phidlam2UTM(Angle phi, Angle dlam)
                        throws Msg {
                return phidlam2GK(phi,dlam).times(0.9996);
        }

/** test program */
  public static void demo() {
                System.out.println("**************************");
                System.out.println("*** Test of Class 'MB' ***");
                System.out.println("**************************");
                System.out.println();

                try {
                Ellipsoid ell = Ellipsoid.WGS84;
    MB b = new MB(ell);
    System.out.println(ell);
    System.out.println(b);
    System.out.println("circumference of meridian ellipse="+
      (b.circumference()/1000.0)+" km");
    //double phi = (47. + 7./60. + 10.19550/3600.)/RHODEG;
    Angle phi = new Angle(47,7,10.19550);
    phi = new Angle(52.,Angle.DEG);
    double G = b.arc(phi);
    System.out.println("arc("+phi.deg()+" deg)="+G+" m");
    Angle phinew = b.lat(G);
    System.out.print("lat("+G+" m)="+phinew.deg()+" deg");
    System.out.println("; (error= "+
        String.format("%.1E",(phi.deg()-phinew.deg())*3600.0)+" arcsec)");
                System.out.println("\nconvergency behaviour");
    for(int n=1;n<=6;n++) {
                        double Gn = new MB(ell,n).arc(phi);
                        System.out.println("N="+n+": arc= "+Gn+" m, error="+
                         (G-Gn)+" m");
                }

                System.out.println("\ncomplex calls to arc/lat");
                Complex z = new Complex(0.9,0.05);
                Complex w = b.arc(z);
                System.out.println("\narc("+z+" rad)="+w+" m");

                phi = new Angle(47,4,30.312);
                Angle dlam = new Angle(0,41,41.768);
                Complex phidlam = new Complex(phi.rad(),dlam.rad());
                System.out.println("\ncomplex "+b+"\nphi="+phi+", delta lambda="+dlam);
                Complex UTM = b.phidlam2UTM(phi,dlam);
                System.out.println("UTM: northing="+UTM.re()+" m, easting="+UTM.im()+" m");
                }
                catch(Msg m) { System.err.println(m); }
                catch(Exception ex) { System.err.println("oops, unknown exception"); }
/*
**************************
*** Test of Class 'MB' ***
**************************

meridian arcs on Bessel, order N=8
circumference of meridian ellipse=40003.42305774125 km
arc(47.11949875 deg)=5220000.551362536 m
lat(5220000.551362536 m)=47.11949874999999 deg; (error= 2.6E-11 arcsec)

convergency behaviour
N=1: arc= 5220020.237369137 m, error=-19.68600660096854 m
N=2: arc= 5220000.608289172 m, error=-0.05692663602530956 m
N=3: arc= 5220000.551573537 m, error=-2.1100137382745743E-4 m
N=4: arc= 5220000.55136343 m, error=-8.940696716308594E-7 m
N=5: arc= 5220000.551362541 m, error=-4.6566128730773926E-9 m
N=6: arc= 5220000.551362536 m, error=0.0 m
*/
  } // end demo()

// public static void main(String [] args) { MB.demo(); }

} // end class MB