/* ----------------------------------------
 * Origin  : GFI INFORMATIQUE 
 * Project : CDPP Propagation Tool 
 * ----------------------------------------
 * History 
 * Version : 01/05/2015 : SCA : V2.1
 * Version : 01/09/2016 : SCA : V2.2 (PTSW H2020)
 *  Computation only in real time.  
 * ----------------------------------------
 */
package fr.tlse.irap.tool.radial;

import fr.tlse.irap.utils.CoordConverter;

public class DragModel {

	private static double rSunM = CoordConverter.rSUN * 1e3; // solar radius in meter
	private static int cd = 1; 					// Dimensionless drag coefficient 
	private static double mi = 1.67 * 1e-27; 	// Mass of proton in kg
	private static double mCME = 1e12; 			// Typical mass of CME

	private int res; 	// time resolution (expressed in seconds) 
	private int dt;	 	// time duration (expressed in seconds)
	private int steps;  // nb steps 
	private long ti;	// initial time of deceleration phase
	private long vi;	// initial constant speed solar wind
	private long vw;
	private double riKm;
	private double gamma;

	public DragModel() {
	}

	/**
	 * @param duration in days
	 * @param timeRes in minutes
	 */
	public DragModel(int duration, int timeRes) {
		this.dt = 3600*24*duration;
		this.res = 60*timeRes;
		this.steps = dt/res;
	}

	/**
	 * function[time,rt,vt]=Dragmodel(speed, extentCME, innerBoundary, ambiantSpeed, dragParam)
	 * <p>Output:
	 * <ul><li>time[steps] :<> in seconds from launch time at the Sun (access via getTime(step))
	 * <li>rt[steps] in kms (access via getRadius(step))
	 * <li>vt[steps] speed of CME in km/s along trajectory (access via getSpeed(step))
	 * <li>gamma : computed drag param (in km^-1)</ul>
	 * @param speed : speed initial speed (km/s) of CME at the Sun
	 * @param extentCME : extent of CME in degrees (Extent of CME in current propagation tool)
	 * @param innerBoundary : expressed in solar radii, (defines the distance from which the Sun accelerates or deccelerates)
	 * @param ambiantSpeed : ambiant solar wind speed (km/s)
	 * @param dragParam : (in 10^-7 km^-1) defines efficiency of drag<dt>
	 */
	public void compute(double speed, double extentCME, double innerBoundary, double ambiantSpeed, double dragParam) {
		
		long ri = (long)innerBoundary;
		double phiCME = Math.toRadians(extentCME);
		riKm = CoordConverter.rsunToKm(ri);
		vw = (long)ambiantSpeed;
		
		// Computed initial time of deceleration phase 
		vi = (long)speed;
		ti = (long)(riKm/vi);

		// Define solar wind model (Leblanc et al. 1998) for nw expressed in m^-3:
		double nw = 1e11*((2.8*Math.pow(ri,-2))+(35*Math.pow(ri,-4))+(680*Math.pow(ri,-6)));
		
		// Calculate the area of the CME:
		double aCME = Math.pow((rSunM*ri*phiCME)/2, 2); // expressed in m^2
		if (dragParam == 0) {
			// Calculate gamma (drag coefficient):
			gamma = (nw*cd*mi*aCME*1e3)/mCME; // expressed in km^-1
		}
		else {
		    gamma = dragParam*1e-7;
		}
	}

	/**
	 * Compute speed from time (in km/s)
	 * @param tDecc
	 * @return
	 */
	private double computeSpeed(long tDecc) {
		if (vi >= vw) {
			return ((vi-vw)/(1+(gamma*(vi-vw)*tDecc)))+vw; // expressed in km/s
		}
		else {
			return ((vi-vw)/(1-(gamma*(vi-vw)*tDecc)))+vw; // expressed in km/s
		}
	}

	/**
	 * Compute distance from time (in km)
	 * @param tDecc : descelerates time
	 * @return
	 */
	private double computeRadius(long tDecc) {
	 	if (vi >= vw) {
	 		return (1./gamma)*Math.log(1+(gamma*(vi-vw)*tDecc))+(vw*tDecc)+riKm; 
		}
		else {
			return -(1./gamma)*Math.log(1-(gamma*(vi-vw)*tDecc))+(vw*tDecc)+riKm; // distances expressed in km, time in seconds
		}
	}

	/**
	 * Get drag time near from radius 
	 * @param radius
	 * @return -1 if not found
	 */
	public long getDragTime(double radius) {
		int step = getStep(radius);
		if (step == 0) {
			// compute propagating time without drag
			return (long)(radius/vi);
		}
		else if (step > 0) {
			// get propagating drag time
			return getTime(step);
		}
		else {
			// not found
			return -1;
		}
	}

	/**
	 * Get radius from start time
	 * @param time in seconds 
	 * @return
	 */
	public double getRadius(long time) {
		double radius=0;
		int step = getStep(time);
		if (step == 0) {
			radius = time*vi;
		}
		else if (step > 0) {
			radius = getRadius(step);
		}
		else {
			radius = computeRadius(time - ti);
		}
		return radius;
	}

	/**
	 * Get speed from start time
	 * @param time in seconds 
	 * @return
	 */
	public double getSpeed(long time) {
		double speed=0;
		int step = getStep(time);
		if (step == 0) {
			speed = vi;
		}
		else if (step > 0) {
			speed = getSpeed(step);
		}
		else {
			speed = computeSpeed(time - ti);
		}
		return speed;
	}

	/**
	 * Get real time in seconds at step
	 * initial time (for constant phase) + computed time (for variable phase)
	 * @param step
	 * @return
	 */
	public long getTime (int step) {
		return ti + getTimeDecc(step);
	}

	/**
	 * Get step from real time
	 * @param time
	 * @return -1 if not found
	 */
	public int getStep(long time) {
		int step = (int)((time - ti)/res);
		if (step < 0) {
			return 0;
		}
		else if (step < steps) {
			return step;
		}
		else {
			return -1;
		}
	}

	/**
	 * Get computed Solar Wind Speed at step 
	 * @param step
	 * @return vt[step]
	 */
	public double getSpeed (int step) {
		return computeSpeed(getTimeDecc(step));
	}

	/**
	 * Get computed Radius at step 
	 * @param step
	 * @return rt[step];
	 */
	public double getRadius (int step) {
		return computeRadius(getTimeDecc(step));
	}

	/**
	 * Get computed Drag Param (10^-7 km-1)
	 * @return
	 */
	public double getDragParam() {
		return gamma*1e7;
	}

	/**
	 * Get time of decceleration in seconds at step
	 * @param step
	 * @return
	 */
	private long getTimeDecc (int step) {
		return res*step;
	}

	/**
	 * Get time duration in seconds
	 * @return
	 */
	public long getDuration() {
		return dt;
	}

	/**
	 * Get time resolution in seconds
	 * @return
	 */
	public long getResolution() {
		return res;
	}

	/**
	 * Get drag time near from radius over the max step 
	 * @param radius
	 * @return -1 if not found
	 */
	public long getDragTimeOver(double radius) {
		int step = getStepOver(radius);
		if (step > 0) {
			// get propagating drag time
			return getTime(step);
		}
		else {
			// not found
			return -1;
		}
	}

	/**
	 * Get Step near Radial distance 
	 * @param radius
	 * @return -1 if not found
	 */
	public int getStep(double radius) {
		double rti0=0;
		double rti1=0;
		for (int i=0; i<steps; i++) {
			rti0 = rti1;
			rti1 = computeRadius(getTimeDecc(i));
			if (rti1 >= radius) {
				if (i > 0) {
					if (radius < Math.abs((rti1-rti0)/2)) {
						return i-1;
					}
					else {
						return i;
					}
				}
				else {
					return i;
				}
			}
		}
		return -1;
	}

	/**
	 * Get Step near Radial distance over the max step 
	 * @param radius
	 * @return -1 if not found
	 */
	public int getStepOver(double radius) {
		int i=steps;
		double rti=0;
		boolean found=false;
		while (!found) {
			rti = computeRadius(getTimeDecc(i));
			if (rti >= radius) {
				found = true;
			}
			else if (i > 30000) {
				// limit the search to avoid infinity loop 
				found = true;
			}
			i++;
		}
		return i;
	}

	/*
	 * Accessors needed for WS data persistence
	 */
	public int getRes() {
		return res;
	}
	public void setRes(int res) {
		this.res = res;
	}
	public int getDt() {
		return dt;
	}
	public void setDt(int dt) {
		this.dt = dt;
	}
	public int getSteps() {
		return steps;
	}
	public void setSteps(int steps) {
		this.steps = steps;
	}
	public long getTi() {
		return ti;
	}
	public void setTi(long ti) {
		this.ti = ti;
	}
	public long getVi() {
		return vi;
	}
	public void setVi(long vi) {
		this.vi = vi;
	}
	public long getVw() {
		return vw;
	}
	public void setVw(long vw) {
		this.vw = vw;
	}
	public double getRiKm() {
		return riKm;
	}
	public void setRiKm(double riKm) {
		this.riKm = riKm;
	}
	public double getGamma() {
		return gamma;
	}
	public void setGamma(double gamma) {
		this.gamma = gamma;
	}
}
