path: root/org/madore/ephem/Comput.java

package org.madore.ephem;

import java.util.List;
import java.util.HashMap;
import java.util.Arrays;

public final class Comput {

    public static final double degree = Math.PI/180;
    public static final double arcsecond = degree/3600;

    public static final double ipow(double x, int n) {
	if ( n == 0 )
	    return 1;
	else if ( n == 1 )
	    return x;
	else if ( n < 0 ) {
	    x = 1./x;
	    n = -n;
	}
	double y = x;
	int k = 1;
	double z = 1;
	int m = 0;
	while ( true ) {
	    if ( ( n & k ) != 0 ) {
		z *= y;
		m += k;
	    }
	    if ( m < n ) {
		y *= y;
		k *= 2;
	    } else
		return z;
	}
    }

    private static Double[] boxDoubleArray(double[] c) {
	Double[] boxed = new Double[c.length];
	for ( int i=0 ; i<boxed.length ; i++ )
	    boxed[i] = c[i];
	return boxed;
    }

    public static interface Function {
	public double v(double t);
    }

    public static final class Constant implements Function {
	double val;
	public Constant(double val) {
	    this.val = val;
	}
	public double v(double t) {
	    return val;
	}
    }

    public static final class Identity implements Function {
	public double v(double t) {
	    return t;
	}
    }

    public static class SumFunction implements Function {
	List<? extends Function> tms;
	public SumFunction(List<? extends Function> tms) {
	    this.tms = tms;
	}
	public SumFunction(Function[] tms) {
	    this(Arrays.asList(tms));
	}
	public double v(double t) {
	    double z = 0;
	    for ( int i=0 ; i<tms.size() ; i++ )
		z += tms.get(i).v(t);
	    return z;
	}
    }

    public static class LinCombFunction implements Function {
	List<Double> c;
	List<Function> tms;
	public LinCombFunction(List<Double> c, List<Function> tms) {
	    if ( c.size() != tms.size() )
		throw new IllegalArgumentException();
	    this.c = c;
	    this.tms = tms;
	}
	public LinCombFunction(List<Double> c, Function[] tms) {
	    this(c, Arrays.asList(tms));
	}
	public LinCombFunction(Double[] c, List<Function> tms) {
	    this(Arrays.asList(c), tms);
	}
	public LinCombFunction(Double[] c, Function[] tms) {
	    this(Arrays.asList(c), Arrays.asList(tms));
	}
	public LinCombFunction(double[] c, List<Function> tms) {
	    this(boxDoubleArray(c), tms);
	}
	public LinCombFunction(double[] c, Function[] tms) {
	    this(boxDoubleArray(c), tms);
	}
	public double v(double t) {
	    double z = 0;
	    for ( int i=0 ; i<tms.size() ; i++ )
		z += c.get(i) * tms.get(i).v(t);
	    return z;
	}
    }

    public static final class TimesPower implements Function {
	int deg;
	Function fn;
	public TimesPower(int deg, Function fn) {
	    this.deg = deg;
	    this.fn = fn;
	}
	public double v(double t) {
	    return ipow(t,deg) * fn.v(t);
	}
    }

    public static final class Polynomial implements Function {
	List<Double> c;
	public Polynomial(List<Double> c) {
	    this.c = c;
	}
	public Polynomial(Double[] c) {
	    this(Arrays.asList(c));
	}
	public Polynomial(double... c) {
	    this(boxDoubleArray(c));
	}
	public double v(double t) {
	    double z = 0;
	    double y = 1;
	    for ( int i=0 ; i<c.size() ; i++ ) {
		z += c.get(i) * y;
		y *= t;
	    }
	    return z;
	}
    }

    public static class Poisson2Term implements Function {
	// Returns c*t^d*cos(phi+sum(om[i]*fn[i](t)))
	int deg;
	double c;
	double phi;
	double[] om;
	Function[] fn;
	Poisson2Term(int deg, double c, double phi, double[] om, Function[] fn) {
	    if ( om.length != fn.length )
		throw new IllegalArgumentException();
	    this.c = c;
	    this.deg = deg;
	    this.phi = phi;
	    this.om = om;
	    this.fn = fn;
	}
	public static Poisson2Term cphi(int deg, double c, double phi, double[] om, Function[] fn) {
	    return new Poisson2Term(deg, c, phi, om, fn);
	}
	public static Poisson2Term ab(int deg, double a, double b, double om[], Function[] fn) {
	    double c = Math.sqrt(a*a + b*b);
	    double phi = Math.atan2(-b, a);
	    return new Poisson2Term(deg, c, phi, om, fn);
	}
	public double v(double t) {
	    double arg = 0;
	    for ( int i=0 ; i<om.length ; i++ )
		arg += om[i]*fn[i].v(t);
	    return c * ipow(t, deg) * Math.cos(arg + phi);
	}
    }

    public static class SumPoisson2Terms extends SumFunction
	implements Function {
	List<? extends Poisson2Term> tms;
	public SumPoisson2Terms(List<? extends Poisson2Term> tms) {
	    super(tms);
	    this.tms = tms;
	}
	public SumPoisson2Terms(Poisson2Term[] tms) {
	    this(Arrays.asList(tms));
	}
	public double v(double t) {
	    double z = 0;
	    final int npow = 6;
	    double[] tpow = new double[npow];
	    for ( int d=0 ; d<npow ; d++ )
		tpow[d] = ipow(t,d);
	    HashMap<Function,Double> fncache = new HashMap<Function,Double>();
	    for ( int i=0 ; i<tms.size() ; i++ ) {
		Poisson2Term tm = tms.get(i);
		int deg = tm.deg;
		double tp = deg<npow ? tpow[deg] : ipow(t,deg);
		double arg = 0;
		for ( int j=0 ; j<tm.om.length ; j++ ) {
		    if ( ! fncache.containsKey(tm.fn[j]) ) {
			double v = tm.fn[j].v(t);
			fncache.put(tm.fn[j], v);
		    }
		    arg += tm.om[j]*fncache.get(tm.fn[j]);
		}
		z += tm.c * tp * Math.cos(arg + tm.phi);
	    }
	    return z;
	}
    }

    public static final class PoissonTerm extends Poisson2Term implements Function {
	// Returns c*t^d*cos(phi+om*t) = t^d*(a*cos(om*t) + b*sin(om*t))
	double om;
	private static Function identity = new Identity();
	PoissonTerm(int deg, double c, double phi, double om) {
	    super(deg, c, phi, new double[] { om }, new Function[] { identity });
	    this.om = om;
	}
	public static PoissonTerm cphi(int deg, double c, double phi, double om) {
	    return new PoissonTerm(deg, c, phi, om);
	}
	public static PoissonTerm ab(int deg, double a, double b, double om) {
	    double c = Math.sqrt(a*a + b*b);
	    double phi = Math.atan2(-b, a);
	    return new PoissonTerm(deg, c, phi, om);
	}
	public double v(double t) {
	    return c * ipow(t, deg) * Math.cos(om*t + phi);
	}
    }

    public static class SumPoissonTerms extends SumPoisson2Terms
	implements Function {
	List<? extends PoissonTerm> tms;
	public SumPoissonTerms(List<? extends PoissonTerm> tms) {
	    super(tms);
	    this.tms = tms;
	}
	public SumPoissonTerms(PoissonTerm[] tms) {
	    this(Arrays.asList(tms));
	}
	public double v(double t) {
	    double z = 0;
	    final int npow = 6;
	    double[] tpow = new double[npow];
	    for ( int d=0 ; d<npow ; d++ )
		tpow[d] = ipow(t,d);
	    for ( int i=0 ; i<tms.size() ; i++ ) {
		PoissonTerm tm = tms.get(i);
		int deg = tm.deg;
		double tp = deg<npow ? tpow[deg] : ipow(t,deg);
		z += tm.c * tp * Math.cos(tm.om*t + tm.phi);
	    }
	    return z;
	}
    }

}