// ModelEbm.java

// import java.io.*;


/** 
 *  Ebm-Model computing the data.
 * 
 */
class ModelEbm extends Model {

    // EBM-Constants: 
    private final double dSigma = 5.7E-8;
    private final double emissivity = 0.95;

    /* EBM-Parameters:
     */
    // Number of timesteps: 
    public long     nTimesteps = 0 ;
    // Initil temperature: 
    public double   dTempBegin ;
    // timestep dt in days: 
    private double  dDtDays ;
    // heat capacity: 
    private double  dCw;
    // albedo feedback on / off: 
    private boolean bAlbedoFeedback ;
    // average albedo: 
    private double  dAlbedoMean;
    // average transmission: 
    private double  dTransmissMean;
    // short-wave incoming energy flux: 
    private double  dFluxSwIn ;
    // noise on / off:
    private boolean bNoiseOn ;
    // initial integer value of the random generator: 
    private long    lNoiseBegin = 0 ;
    // amplitude of the noise: 
    private double  dNoiseAmplitude ;

    /* Variables used internally for computation: 
     */
    // Timestep in basic unit: seconds: 
    private double  dDt;
    // actual albedo at timestep: 
    private double  dAlbedo;
    // actual transmission at timestep: 
    private double  dTransmiss;
    // short-wave outgoing energy flux:
    private double  dFluxSwOut ;
    // long-wave outgoing energy flux:
    private double  dFluxLwOut ;
    // net flux:
    private double  dFluxnet ;

    public ModelStateEbm state     = new ModelStateEbm() ;    
    public ModelStateEbm stateNext = new ModelStateEbm(); 

    public RandomGaussian randomNoise ; 

    //  Step-counter: needed only for testing: 
    private long i ;

    //PrintWriter printlog = new PrintWriter(System.out, true);



    ModelEbm(){
    }


    /** 
     * Returns the current state
     */
    public ModelState getState (){
        return state ;
    }


    /** 
     * Returns the next state
     */
    public ModelState getStateNext (){
        return stateNext ;
    }


    /** 
     * Returns the number of timesteps
     */
    public long getTimesteps (){
        return nTimesteps;
    }


    /** 
     * Resets the model to initial conditions 
     * (state and random generator varaibles
     */
    public void resetModel() {
        i = 0 ;

	randomNoise = new RandomGaussian(); 
        randomNoise.RandomGaussianInit( lNoiseBegin );

	state.dTemp  = dTempBegin;
    }


    /** 
     * Initializaton of the EbmModel with the parameters
     * and reset. 
     */
    public ModelStateEbm initModel( ParameterSetEbm param )
    {
        // Set parameters:
        nTimesteps =      param.Timesteps.getValue();
        dTempBegin =      param.TempBegin.getValue();
        dDtDays =         param.TimeStepLength.getValue();
        dCw =             param.HeatCapac.getValue();
        bAlbedoFeedback = param.AlbedoFeedback.getValue();
        dAlbedoMean =     param.AlbedoMean.getValue();
        dTransmissMean =  param.TransmissMean.getValue();
        dFluxSwIn =       param.FluxSwIn.getValue();
        bNoiseOn =        param.NoiseOn.getValue();
        lNoiseBegin =     param.NoiseBegin.getValue() * (-1);
        dNoiseAmplitude = param.NoiseAmplitude.getValue();   

        // Internal constants and transformations: 
	dDt = 60*60*24 * (double)dDtDays ;

	//        randomNoise = new RandomGaussian(); 

        // Reset state and internal varaiables (here: noise)
        resetModel(); 

	return state ;
    }



    /** 
     * Perform one timstep, i.e. compute the state of next timestep
     * (i+1) from the current state (i), 
     * and set the new current state to the values of i+1. 
     */
    public void stepModel()
    {
	// Returns state, i.e. temperature, for step+1

	// Albedo: with or without feedback
	if ( bAlbedoFeedback == true ){
            dAlbedo = dAlbedoMean - dAlbedoMean * 
	        0.025 * tanh( 1.548 * ( state.dTemp - 288.0 ) ) ;
        }
	else{
	    dAlbedo = dAlbedoMean;
	};
        

        // Transmission: with or without noise:
        double dTransmissNoise; 
	if ( bNoiseOn == true ){                

            state.dNoise = randomNoise.NextRandomGaussian();
            dTransmissNoise = dNoiseAmplitude * state.dNoise;
               ;

	    dTransmiss = dTransmissMean + dTransmissMean * dTransmissNoise;

	}
	else{
	    dTransmiss = dTransmissMean;
            state.dNoise = 0.0;
	};

	// Fluxes: 
	dFluxSwOut = dAlbedo * dFluxSwIn;
	dFluxLwOut = emissivity * dSigma * 
	    state.dTemp * state.dTemp * state.dTemp * state.dTemp ;

        dFluxnet = dFluxSwIn - dFluxSwOut - dFluxLwOut * dTransmiss ;
		
	// Temperature for next time-step
	stateNext.dTemp = state.dTemp + dDt * (1/dCw) * dFluxnet ;

        state.setEqual( stateNext );

    }



	public double tanh( double dArgument ){
	// Returns the hyperbolic tangent of dArgument
        if ( dArgument > 25 ) {
		return dArgument / Math.abs( dArgument );
	}
	else{
	    if ( Math.abs(dArgument) < 1.0E-12 ){
	        return dArgument ;
	    }
	    else{
	      	return ( Math.exp(dArgument) - Math.exp(-dArgument) ) /
		   ( Math.exp(dArgument) + Math.exp(-dArgument) ) ;
	        };
	    };
        }


    public String getAspectName ( int iAspectVarIndex ){
        String aspectName ;

        switch ( iAspectVarIndex ) {
            case 0:      // Temperature
	        aspectName = " - Temperature (K)";
                break;        

            default:     // programming error
	        aspectName = "";
                break;        
	};

        return aspectName;          
    }


}