/* ----------------------------------------
 * Origin  : GFI INFORMATIQUE 
 * Project : CDPP Propagation Tool 
 * ----------------------------------------
 * History 
 * Version : 01/05/2015 : SCA : V2.1
 * Version : 07/04/2016 : SCA : V2.2 (PSWS-H2020)
 * => Add Z axis in propagation computation 
 * ----------------------------------------
 */
package fr.tlse.irap.utils;

import fr.tlse.irap.spaceElement.SpaceElement;
import fr.tlse.irap.spaceElement.Sun;
import fr.tlse.irap.tool.sep.Sep.Particle;

public class ToolHelper {

	public final static double rE = 215.095; // Solar radius to 1 AU
	public final static double rC = 1.0; 	 // Radius where the spiral starts (AU)
	
	private final static double protonER = 938.7d;
	private final static double electronER = 0.51d;
	private final static double DTC = 8.3;

	/**
	 * Radial distance away surface of the sun
	 * @param spaceElement
	 * @return
	 */
	public static double computeRadialDistance(SpaceElement e) {
		if (e instanceof Sun) {
			return 0.0;
		} else {
			double radius = Math.sqrt(Math.pow(e.getX(), 2.0) + Math.pow(e.getY(), 2.0) + Math.pow(e.getZ(), 2.0));
			return radius - getRadiusSun();
		}
	}

	/**
	 * Radial distance away surface of the sun projected on XY plane
	 * @param spaceElement
	 * @return
	 */
	public static double computeRadialDistanceXY(SpaceElement e) {
		if (e instanceof Sun) {
			return 0.0;
		} else {
			double radius = Math.sqrt(Math.pow(e.getX(), 2.0) + Math.pow(e.getY(), 2.0));
			return radius - getRadiusSun();
		}
	}

	/**
	 * @param rotationPeriod expressed in days
	 * @return rotation frequency expressed in radians/s
	 */
	public static double getRotationFrequency(double rotationPeriod) {
		return (2 * Math.PI / (rotationPeriod * 3600 * 24));
	}

	/**
	 * Radial distance away surface of the sun
	 * @return
	 */
	public static double getRadiusSun() {
		return CoordConverter.toKm(1.0/ rE);
	}

	/**
	 * @param deltaLong
	 * @param extentCME
	 * @return true if deltaLong is less or equals than extentCME or if deltaLong is greater or equals than 360-extentCME
	 */
	public static boolean isInsideCME(double deltaLong, double extentCME) {
		return (deltaLong <= extentCME || deltaLong >= 360.0 - extentCME);
	}

	/**
	 * Calculate the length of the spiral (see script CODEA of SEPpropagation_new.pdf on 17/01/2013)
	 * @param Speed of solar wind, expressed in km/s-1
	 * @param Radial distance away Sun of probe, expressed in AU
	 * @param Sun rotation period, expressed in day
	 * @return Length of the spiral, expressed in AU
	 */
	public static double computeSpiralLength(double speed, double radialDistance, double rotationPeriod) {

		double vp = speed * 0.124138; 				// Solar wind speed, expressed in solar radii/day
		double timeStep = 0.0005; 					// Time step of the spiral (0.7 minutes), expressed in day
		double omegaRadians = (2 * Math.PI) / rotationPeriod; // Rotation frequency, expressed in radians/day
		Double radiiDistance = radialDistance * rE; // Radial distance away Sun of probe, expressed in solar radii/day
		
		/* Calculate the length of the spiral */
		double length = 0.0d;
		double dr = vp * timeStep;
		double dc;
		double dphi = omegaRadians * timeStep;
		Double radius = rC;

		while (radius.compareTo(radiiDistance) <= 0) {
			radius += dr;
			dc = Math.sqrt(((radius - dr) * (radius - dr)) + (radius * radius)
					- (2 * radius * (radius - dr) * Math.cos(dphi)));
			length += dc;
		}
		return (length / rE);
	}

	/**
	 * Indicates if a probe or planet is inside when spread less than 180°.
	 * 
	 * Test the following relation : r <= rEND <= R where r is associated to the
	 * nearest point of INTERNAL spiral and R is associated to the nearest point
	 * of EXTERNAL spiral and rEND is associated to the given point to check.
	 * All longitudes must be expressed in the same system coordinates (HAE, Carrington, ...)
	 * 
	 * @param speed
	 * @param rotationPeriod
	 * @param fpt1 : internal footpoint longitude 
	 * @param fpt2 : external footpoint longitude
	 * @param longTarget : target longitude 
	 * @param distTarget : target radial distance (AU)
	 * @return
	 */
	public static boolean isInsideSpiral(double speed, double rotationPeriod, double fpt1, double fpt2, double longTarget, double distTarget) {

		double omegaRadians = (2 * Math.PI) / rotationPeriod; // Rotation frequency, expressed in radians/day
		double vp = speed * 0.124138; 				// Solar wind speed, expressed in solar radii/day	
		double coeffA = rC / rE; 					// Coefficient A to turn the spiral, expressed in AU		
		double coeffB = vp / (omegaRadians * rE);	// Coefficient B to determine distance between arms of spiral expressed in AU by radians
		double theta = Math.toRadians(longTarget);
		double theta1 = Math.toRadians(fpt1);
		double theta2 = Math.toRadians(fpt2);
		double rho1 = 0.0d;
		double rho2 = 0.0d;
		boolean isInside=false;
		int i = 0;
		
		while (rho1 < distTarget) {
			// Calculate the point of INTERNAL spiral the nearest to the given point
			theta = theta - (i * 2 * Math.PI);
			rho1 = coeffA + (coeffB * (theta1 - theta));
			if (rho1 <= distTarget) {
				// Calculate the point of EXTERNAL spiral the nearest to the given point
				theta = theta - (i * 2 * Math.PI);
				rho2 = coeffA + (coeffB * (theta2 - theta));
				isInside = (distTarget <= rho2);
				i++;
			}
		}
		return (isInside);
	}
	
	/**
	 * Calculate time of SEP propagation of particle between two points; see script CODEB from SEPpropagation_new.pdf on 17/01/2013
	 * @param particleEnergy : energy of particle
	 * @param particleType : type of particle
	 * @param spiralLength : length of spiral between the two points
	 * @return Time of propagation in minutes
	 */
	public static double computePropagationTime(double particleEnergy, Particle particle, double spiralLength) {
		// Rest energy of particle according to its type
		double restEnergy;
		switch (particle) {
		case proton:
			restEnergy = protonER;
			break;
		case electron:
			restEnergy = electronER;
			break;
		default:
			restEnergy = protonER;
			//System.out.println("unexpected type of particle ==> force to proton");
			break;
		}
		// Total energy of particle
		double totalEnergy = particleEnergy + restEnergy;
		//System.out.println("computePropagationTime with : DTC=" + DTC + " restEnergy=" + restEnergy + " totalEnergy=" + totalEnergy);
		Double deltaMinutes = ((spiralLength * DTC) / Math.sqrt(1 - Math.pow(restEnergy / totalEnergy, 2)));
		if (particleEnergy == 0 || deltaMinutes.isNaN() || deltaMinutes.isInfinite()) {
			deltaMinutes = 0.0d;
			//System.out.println("Propagation time is forced to 0, check particleEnergy=" + particleEnergy);
		}
		return deltaMinutes;
	}

}
