//PlaneWave_State.java
import java.util.*;


public class PlaneWave_State {
    //public static final double epsilon0 = 8.8541878176E-12; //Units: F/m
    // Approximate epsilon0 value consistent with approximate velocity
    public static final double epsilon0 = 8.841941286E-12; //Units: F/m
    public static final double mu0 = 1.25663706144E-6; //Units H/m
    //private static final double c = 1.0/Math.sqrt(epsilon0*mu0);
    private static final double c = 3.0E8; // Approximate value
    private static final double c_light2 = 9.0E16; // Approximate value
    private static final double electron_charge = 1.60217646e-19; // Coulombs
    private static final double electron_mass = 9.10938188e-31; // Kg
    private static final double proton_mass = 1.67262158e-27; // Kg
    private static final double neutron_mass = 1.67492729e-27; // Kg
    
    double phase_velocity, vario, kineticE, velocity, velocity2;
    public double Ex, Ey, Hz, Bz; //Amplitudes of fields
    public double EleFx, EleFy, MFx, MFy;
    
    public double EnergyMax, velocityMax, dt_suggested;
    
    double epsilon_r, mu_r; //relative permeability and permitivity
    public double xpos, ypos, xpos_old, ypos_old, vmax, vmin, tempo,
                  xmin, xmax, ymin, ymax, Deltax, Deltay, Deltavx, Deltavy;
    public double xpos_0, ypos_0, vx_0, vy_0, vx, vy, vx_old, vy_old;
    public double ax, ay, ax_old, ay_old;
    public int SCROLLMAX2 = 20001;
    public int particle, jj;
    
    public double mantissa;
    public int exponent;
    
    public double xnewmin, xnewmax, ynewmin, ynewmax;
    public double charge, mass, charge_value, mass_value;
    
    int NPointsZ, NPointsT, DeltaT, CurrentTime, TotalTime;
    double dt, dx, dy;
    double  Phix, Phiy, wt;
    int current_z;
    public boolean IsExample1, IsExample2, IsExample3, IsGeneral_Input, IsOpen;
    boolean RESET;
    public int SleepTime;
    
       
    public boolean LicenseExpired;
    public int this_month, today_week, this_year, this_hour, this_minute, today_month, 
	       today_year,this_zone, saving_time;
    GregorianCalendar Greg = new GregorianCalendar();

    
    public PlaneWave_State(){    
	//Independent Variables
	epsilon_r = 1.0;
	mu_r = 1.0;
	NPointsZ = 360;
	NPointsT = 360;
	DeltaT = 1; //timestep increment
        dt = 1.0e-9; // seconds (time step)
        IsExample1 = true;
        IsExample2 = false;
        IsExample3 = false;
        IsGeneral_Input = false;
        IsOpen = false;
	
        tempo = 0.0; // actual cumulative time in seconds
        particle = 1; // electron
        jj = 0;
        
        charge = -1.0;
        charge_value = charge * electron_charge;
        
        mass = 1.0;
        mass_value = mass * electron_mass;
        
        CurrentTime = 0;
        TotalTime = 0;
        SleepTime = 50;
        
	Ex = 1.0; // V/m
        Ey = 1.0; // V/m
	Hz = 1.0; // A/m       Ex/Math.sqrt(epsilon0*mu0);
	Bz = Hz*mu0;
        
        vario = charge_value/mass_value*dt;
        
        Phix = 0.0;
	current_z = 0;
	RESET = false;
	ignition(); 
        xpos = 0.0;
        ypos = 0.0;
        xpos_0 = 0.0;
        ypos_0 = 0.0;
        xpos_old = 0.0;
        ypos_old = 0.0;
        
        vx_0 = 1.0;
        vy_0 = 1.0;
        vx = vx_0*1.0e5;
        vy = vy_0*1.0e5;
        vx_old = vx;
        vy_old = vy;
        
        ax_old = charge_value/mass_value*(Ex + vy_old * Bz);
        ay_old = charge_value/mass_value*(Ey - vx_old * Bz);
        ay = ay_old;
        ax = ax_old;
        
        velocity2 = vx*vx+vy*vy;
        velocity = Math.sqrt(velocity2);
        kineticE = 0.5*mass_value*velocity2;
        
        ax = 0.0;
        ay = 0.0;
        
        //xmin = -1.0;
        //xmax = 1.0;
        //ymin = -1.0;
        //ymax = 1.0;
        
        xmin = -0.1;
        xmax = 0.1;
        ymin = -0.1;
        ymax = 0.1;
        
        
        vmin = -10.0;
        vmax = 10.0; 
                
        Deltax = 2.0*xmax/(double)(SCROLLMAX2-1);
        Deltavx = 2.0*vmax/(double)(SCROLLMAX2-1);
        
        //Deltay = 2.0*ymax/(double)(SCROLLMAX2-1);
        
	this_month = Greg.get(Calendar.MONTH);
	today_week = Greg.get(Calendar.DAY_OF_WEEK);
	this_year = Greg.get(Calendar.YEAR);
	this_hour = Greg.get(Calendar.HOUR_OF_DAY);
	this_minute = Greg.get(Calendar.MINUTE);
	today_month = Greg.get(Calendar.DAY_OF_MONTH);
	today_year = Greg.get(Calendar.DAY_OF_YEAR);
	this_zone = Greg.get(Calendar.ZONE_OFFSET);
	saving_time = Greg.get(Calendar.DST_OFFSET);
    }
    
    public synchronized void ignition(){
	//Dependent Variables
        Deltax = 2.0*xmax/(double)(SCROLLMAX2-1);
        Deltavx = 2.0*vmax/(double)(SCROLLMAX2-1);
        
        Bz = Hz*mu0;
        vx = vx_0*1.0e5;
        vy = vy_0*1.0e5;
        vx_old = vx;
        vy_old = vy;
        
        velocity2 = vx*vx+vy*vy;
        velocity = Math.sqrt(velocity2);
        kineticE = 0.5*mass_value*velocity2;
        
        if(particle == 1){// electron
            charge = -1.0;
            charge_value = charge * electron_charge;
            mass_value = electron_mass;
            mass_value = electron_mass/Math.sqrt(1-velocity2/c_light2);
        }
        else if(particle == 2){// proton
            charge = 1.0;
            charge_value = charge * electron_charge;
            mass_value = proton_mass;
            mass_value = proton_mass/Math.sqrt(1-velocity2/c_light2);
        }
        else if(particle == 3){// neutron
            charge = 0.0;
            charge_value = charge * electron_charge;
            mass_value = neutron_mass;
            mass_value = neutron_mass/Math.sqrt(1-velocity2/c_light2);
        }
       
        ax = charge_value/mass_value*(Ex + vy_old * Bz);
        ay = charge_value/mass_value*(Ey - vx_old * Bz);
        
        EleFx = charge_value/mass_value*Ex;
        EleFy = charge_value/mass_value*Ey;
        MFx = charge_value/mass_value*(vy_old * Bz);
        MFy = charge_value/mass_value*( - vx_old * Bz);
        
        //System.out.println(" Mx =  "+MFx+" My =  "+MFy);
        
        double testE, testM;
        
        testE = Math.sqrt(EleFx*EleFx + EleFy*EleFy);
        testM = Math.sqrt(MFx*MFx + MFy*MFy);
        
        vario = charge_value/mass_value*dt;
        
        //EnergyMax = 0.5 * mass_value * velocity2 + Math.abs(charge_value) * 2.0 * xmax * Math.sqrt(2) * (Ex*Ex + Ey*Ey);
        EnergyMax = 0.5 * mass_value * (1.4e5*1.4e5) + Math.abs(charge_value) * 2.0 * xmax * Math.sqrt(2) * (Ex*Ex + Ey*Ey);
        velocityMax = Math.sqrt(2.0*EnergyMax/mass_value);
        double multiplier = 1.0;
        double dtE = 1.0, dtM = 1.0;
        
            dtE = 0.001/velocityMax;
            if(Math.abs(Hz) <= 10.0){multiplier = 0.01;}
            else{
                multiplier = 1.0/(100.0 + (Math.abs(Hz)-10.0)/(990.0/40.0));  
            }
            dtM = 0.5/Math.abs(velocityMax*Hz)*multiplier;
            
        if(testM == 0.0){
            dt_suggested = dtE;
        }
        else if(testE == 0.0){
            
            dt_suggested = dtM;
        }
        else{
            dt_suggested = Math.min(dtE, dtM);
        }
        
        if(particle == 2){
            dt_suggested = 1.0e3*dt_suggested;
        }
            
            mantissa = 0.0;
            exponent = 0;
            double test;
            test = dt_suggested;
            
            if(test >= 1.0e-1){
                mantissa = 1.0;
                exponent = -1;
            }
            else if(test < 1.0e-1 && test >= 1.0e2){
                mantissa = MaestroA.rounder(test*1.0e2,3);
                exponent = -2;
            }
            else if(test < 1.0e-2 && test >= 1.0e-3){
                mantissa = MaestroA.rounder(test*1.0e3,3);
                exponent = -3;
            }
            else if(test < 1.0e-3 && test >= 1.0e-4){
                mantissa = MaestroA.rounder(test*1.0e4,3);
                exponent = -4;
            }
            else if(test < 1.0e-4 && test >= 1.0e-5){
                mantissa = MaestroA.rounder(test*1.0e5,3);
                exponent = -5;
            }
            else if(test < 1.0e-5 && test >= 1.0e-6){
                mantissa = MaestroA.rounder(test*1.0e6,3);
                exponent = -6;
            }
            else if(test < 1.0e-6 && test >= 1.0e-7){
                mantissa = MaestroA.rounder(test*1.0e7,3);
                exponent = -7;
            }
            else if(test < 1.0e-7 && test >= 1.0e-8){
                mantissa = MaestroA.rounder(test*1.0e8,3);
                exponent = -8;
            }
            else if(test < 1.0e-8 && test >= 1.0e-9){
                mantissa = MaestroA.rounder(test*1.0e9,3);
                exponent = -9;
            }
            else if(test < 1.0e-9 && test >= 1.0e-10){
                mantissa = MaestroA.rounder(test*1.0e10,3);
                exponent = -10;
            }
            else if(test < 1.0e-10 && test >= 1.0e-11){
                mantissa = MaestroA.rounder(test*1.0e11,3);
                exponent = -11;
            }
            else if(test < 1.0e-11 && test >= 1.0e-12){
                mantissa = MaestroA.rounder(test*1.0e12,3);
                exponent = -12;
            }
            else if(test < 1.0e-12 && test >= 1.0e-13){
                mantissa = MaestroA.rounder(test*1.0e13,3);
                exponent = -13;
            }
            else if(test < 1.0e-13 && test >= 1.0e-14){
                mantissa = MaestroA.rounder(test*1.0e14,3);
                exponent = -14;
            }
            else if(test < 1.0e-14 && test >= 1.0e-15){
                mantissa = MaestroA.rounder(test*1.0e15,3);
                exponent = -15;
            }
            else if(test < 1.0e-15 && test >= 1.0e-16){
                mantissa = MaestroA.rounder(test*1.0e16,3);
                exponent = -16;
            }
            else if(test < 1.0e-16 && test >= 1.0e-17){
                mantissa = MaestroA.rounder(test*1.0e17,3);
                exponent = -17;
            }
            else if(test < 1.0e-17 && test >= 1.0e-18){
                mantissa = MaestroA.rounder(test*1.0e18,3);
                exponent = -18;
            }
            else if(test < 1.0e-18 && test >= 1.0e-19){
                mantissa = MaestroA.rounder(test*1.0e19,3);
                exponent = -19;
            }
            else if(test < 1.0e-19 && test >= 1.0e-20){
                mantissa = MaestroA.rounder(test*1.0e20,3);
                exponent = -20;
            }
            else if(test < 1.0e-20){
                mantissa = 1.0;
                exponent = -20;
            }
            
        //System.out.println("EnergyMax = "+EnergyMax+"   velMax = "+velocityMax+"  mult = "+multiplier);
    }
    
    public synchronized void increment(){
	double testx, testy, dxout, dyout, dtout, vxout, vyout;
        tempo = tempo + dt;
        jj+=1;
        
        velocity2 = vx*vx+vy*vy;
        if(velocity2 > 1.0E16){ //relativistic correction for mass in case particle gets too fast
            if(particle == 1){// electron
                mass_value = electron_mass/Math.sqrt(1-velocity2/c_light2);
                vario = charge_value/mass_value*dt;
            }
            else if(particle == 2){// proton
                mass_value = proton_mass/Math.sqrt(1-velocity2/c_light2);
                vario = charge_value/mass_value*dt;
            }
            else if(particle == 3){// neutron
                mass_value = neutron_mass/Math.sqrt(1-velocity2/c_light2);
                vario = 0.0;
            }
        }
        
        dxout = 0.0;
        dyout = 0.0;
        vxout = 0.0;
        vyout = 0.0;
        dtout = dt;
        
        xpos_old = xpos; 
        ypos_old = ypos;
        vx_old = vx;  
        vy_old = vy;
        
        // definition of accelerations - used for Runge-Kutta;
        ax = charge_value/mass_value*(Ex + vy_old * Bz);
        ay = charge_value/mass_value*(Ey - vx_old * Bz);
        
        // These are used in PlaneWaveCanvas to plot forces at the particle
        EleFx = charge_value/mass_value*Ex;
        EleFy = charge_value/mass_value*Ey;
        MFx = charge_value/mass_value*(vy_old * Bz);
        MFy = charge_value/mass_value*( - vx_old * Bz);
        
        //simple integration
        //vx = vx_old + vario*(Ex + vy_old * Bz);
        //vy = vy_old + vario*(Ey - vx_old * Bz);
        //dx = dt * vx; dy = dt * vy;
        
        //Runge-Kutta (working)
        double dx2, vx2, dy2, vy2;
        double K1, K2, K3, K4;
        K1 = dt*vx;
        K2 = dt*vx + 0.5*K1;
        K3 = dt*vx + 0.5*K2;
        K4 = dt*vx + K3;
        dx = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);
        
        K1 = dt*ax;
        K2 = dt*ax + 0.5*K1;
        K3 = dt*ax + 0.5*K2;
        K4 = dt*ax + K3;
        vx = vx_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);
        
        K1 = dt*vy;
        K2 = dt*vy + 0.5*K1;
        K3 = dt*vy + 0.5*K2;
        K4 = dt*vy + K3;
        dy = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);
        
        K1 = dt*ay;
        K2 = dt*ay + 0.5*K1;
        K3 = dt*ay + 0.5*K2;
        K4 = dt*ay + K3;
        vy = vy_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);
        
        // Corrector pass
        // note average of velocity to improve conservation of energy in magnetic term
        // I invented this!!!
        
        ax = charge_value/mass_value*(Ex + (vy_old+vy)*0.5 * Bz);
        ay = charge_value/mass_value*(Ey - (vx_old+vx)*0.5 * Bz);
        
        double vx1;
        double rate = 0.35;
        vx1 = (rate*vx_old+(1.0-rate)*vx);
        K1 = dt*vx1;
        K2 = dt*vx1 + 0.5*K1;
        K3 = dt*vx1 + 0.5*K2;
        K4 = dt*vx1 + K3;
        dx = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);
        
        K1 = dt*ax;
        K2 = dt*ax + 0.5*K1;
        K3 = dt*ax + 0.5*K2;
        K4 = dt*ax + K3;
        vx = vx_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);
        
        double vy1;
        vy1 = (rate*vy_old+(1.0-rate)*vy);
        
        K1 = dt*vy1;
        K2 = dt*vy1 + 0.5*K1;
        K3 = dt*vy1 + 0.5*K2;
        K4 = dt*vy1 + K3;
        dy = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);
        
        K1 = dt*ay;
        K2 = dt*ay + 0.5*K1;
        K3 = dt*ay + 0.5*K2;
        K4 = dt*ay + K3;
        vy = vy_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);
        
        //----------------------------------------------------------------------
        //System.out.println(jj+"   "+dx+"   "+dx2+"   "+vx+"   "+vx2);
        
        velocity2 = vx*vx+vy*vy;
        velocity = Math.sqrt(velocity2);
        kineticE = 0.5*mass_value*velocity2;
        
        testx = xpos + dx;
        testy = ypos + dy;
        
        if(IsOpen){ // Open Boundaries
                xpos = testx;  ypos = testy;
        }
        else{ // Closed Boundaries
            // particle stays inbounds
            if(Math.abs(testx) < xmax && Math.abs(testy) < ymax){ 
                xpos = testx;  ypos = testy;
            }
            else{ // particle bounces
                if(Math.abs(testx) < xmax && (testy) > ymax){
                    if(ax == 0.0 && ay == 0.0){ // no acceleration, simple reflection
                        xpos = testx;
                        ypos = ymax - (testy-ymax); //ypos - dy;
                        vy=-vy;
                    }
                    else{// use simple Euler to determine time of exit
                        //xpos = testx;
                        dyout = dy - (testy - ymax);
                        vyout = vy_old + vario*(Ey - vx_old * Bz);
                        dtout = dyout/vyout; // time  boundary
                        //vy = -vyout + vario*(Ey - vx_old * Bz);
                        //dy = (dt-dtout) * vy;
                        //ypos = ymax -dy;
                        
                        ax = charge_value/mass_value*(Ex + (vy_old+vy)*0.5 * Bz);
                        ay = charge_value/mass_value*(Ey - (vx_old+vx)*0.5 * Bz);

                        K1 = dtout*vx;
                        K2 = dtout*vx + 0.5*K1;
                        K3 = dtout*vx + 0.5*K2;
                        K4 = dtout*vx + K3;
                        dx = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ax;
                        K2 = dtout*ax + 0.5*K1;
                        K3 = dtout*ax + 0.5*K2;
                        K4 = dtout*ax + K3;
                        vx = vx_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*vy;
                        K2 = dtout*vy + 0.5*K1;
                        K3 = dtout*vy + 0.5*K2;
                        K4 = dtout*vy + K3;
                        dy = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ay;
                        K2 = dtout*ay + 0.5*K1;
                        K3 = dtout*ay + 0.5*K2;
                        K4 = dtout*ay + K3;
                        vy = vy_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        xpos = xpos + dx;
                        ypos = ymax;
                        vy = -vy;
                        tempo = tempo - (dt-dtout);
                    }
                }
                else if(Math.abs(testx) < xmax && (testy) < ymin){
                    if(ax == 0.0 && ay == 0.0){ // no acceleration, simple reflection
                        xpos = testx;
                        ypos = ymin + Math.abs(testy-ymin); //ypos + dy;
                        vy=-vy;
                    }
                    else{
                        //xpos = testx;
                        dyout = dy - (testy - ymin );
                        vyout = vy_old + vario*(Ey - vx_old * Bz);
                        dtout = dyout/vyout; // time  boundary
                        //vy = -vyout + vario*(Ey - vx_old * Bz);
                        //dy = (dt-dtout) * vy;
                        //ypos = ymin - dy;
                        
                        ax = charge_value/mass_value*(Ex + (vy_old+vy)*0.5 * Bz);
                        ay = charge_value/mass_value*(Ey - (vx_old+vx)*0.5 * Bz);

                        K1 = dtout*vx;
                        K2 = dtout*vx + 0.5*K1;
                        K3 = dtout*vx + 0.5*K2;
                        K4 = dtout*vx + K3;
                        dx = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ax;
                        K2 = dtout*ax + 0.5*K1;
                        K3 = dtout*ax + 0.5*K2;
                        K4 = dtout*ax + K3;
                        vx = vx_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*vy;
                        K2 = dtout*vy + 0.5*K1;
                        K3 = dtout*vy + 0.5*K2;
                        K4 = dtout*vy + K3;
                        dy = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ay;
                        K2 = dtout*ay + 0.5*K1;
                        K3 = dtout*ay + 0.5*K2;
                        K4 = dtout*ay + K3;
                        vy = vy_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        xpos = xpos + dx;
                        ypos = ymin;
                        //ypos = ypos_old +dy;
                        tempo = tempo - (dt-dtout);
                        vy = -vy;
                    }
                }
                else if(Math.abs(testy) < ymax && (testx) > xmax){
                    if(ax == 0.0 && ay == 0.0){ // no acceleration, simple reflection
                        xpos = xmax - (testx - xmax); //xpos - dx;
                        ypos = testy;
                        vx=-vx;
                    }
                    else{
                        //ypos = testy;
                        dxout = dx - (testx - xmax);
                        vxout = vx_old + vario*(Ex - vy_old * Bz);
                        dtout = dxout/vxout; // time  boundary
                        //vx = -vxout + vario*(Ex - vy_old * Bz);
                        //dx = (dt-dtout) * vx;
                        //xpos = xmax -dx;
                        
                        ax = charge_value/mass_value*(Ex + (vy_old+vy)*0.5 * Bz);
                        ay = charge_value/mass_value*(Ey - (vx_old+vx)*0.5 * Bz);

                        K1 = dtout*vx;
                        K2 = dtout*vx + 0.5*K1;
                        K3 = dtout*vx + 0.5*K2;
                        K4 = dtout*vx + K3;
                        dx = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ax;
                        K2 = dtout*ax + 0.5*K1;
                        K3 = dtout*ax + 0.5*K2;
                        K4 = dtout*ax + K3;
                        vx = vx_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*vy;
                        K2 = dtout*vy + 0.5*K1;
                        K3 = dtout*vy + 0.5*K2;
                        K4 = dtout*vy + K3;
                        dy = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ay;
                        K2 = dtout*ay + 0.5*K1;
                        K3 = dtout*ay + 0.5*K2;
                        K4 = dtout*ay + K3;
                        vy = vy_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        ypos = ypos + dy;
                        xpos = xmax;
                        tempo = tempo - (dt-dtout);
                        vx = -vx;
                        
                    }
                }
                else if(Math.abs(testy) < ymax && (testx) < xmin){
                    if(ax == 0.0 && ay == 0.0){ // no acceleration, simple reflection
                        xpos = xmin + Math.abs(testx - xmin); //xpos + dx;
                        ypos = testy;
                        vx=-vx;
                    }
                    else{
                        //ypos = testy;
                        dxout = dx - (testx - xmin );
                        vxout = vx_old + vario*(Ex - vy_old * Bz);
                        dtout = dxout/vxout; // time  boundary
                        //vx = -vxout + vario*(Ex - vy_old * Bz);
                        //dx = (dt-dtout) * vx;
                        //xpos = xmin - dx;
                        
                        ax = charge_value/mass_value*(Ex + (vy_old+vy)*0.5 * Bz);
                        ay = charge_value/mass_value*(Ey - (vx_old+vx)*0.5 * Bz);

                        K1 = dtout*vx;
                        K2 = dtout*vx + 0.5*K1;
                        K3 = dtout*vx + 0.5*K2;
                        K4 = dtout*vx + K3;
                        dx = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ax;
                        K2 = dtout*ax + 0.5*K1;
                        K3 = dtout*ax + 0.5*K2;
                        K4 = dtout*ax + K3;
                        vx = vx_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*vy;
                        K2 = dtout*vy + 0.5*K1;
                        K3 = dtout*vy + 0.5*K2;
                        K4 = dtout*vy + K3;
                        dy = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ay;
                        K2 = dtout*ay + 0.5*K1;
                        K3 = dtout*ay + 0.5*K2;
                        K4 = dtout*ay + K3;
                        vy = vy_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        ypos = ypos + dy;
                        xpos = xmin;
                        tempo = tempo - (dt-dtout);
                        vx = -vx;
                    }
                }
                // corners 
                else if((testx) > xmax && (testy) > ymax){
                    if(ax == 0.0 && ay == 0.0){ // no acceleration, simple reflection
                        xpos = xmax - (testx - xmax); //xpos - dx;
                        ypos = ymax - (testy-ymax); //ypos - dy;
                        vx = -vx;
                        vy = -vy;
                    }
                    else{
                        dyout = dy - (testy - ymax);
                        vyout = vy_old + vario*(Ey - vx_old * Bz);
                        dtout = dyout/vyout; // time  boundary
                        
                        ax = charge_value/mass_value*(Ex + (vy_old+vy)*0.5 * Bz);
                        ay = charge_value/mass_value*(Ey - (vx_old+vx)*0.5 * Bz);

                        K1 = dtout*vx;
                        K2 = dtout*vx + 0.5*K1;
                        K3 = dtout*vx + 0.5*K2;
                        K4 = dtout*vx + K3;
                        dx = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ax;
                        K2 = dtout*ax + 0.5*K1;
                        K3 = dtout*ax + 0.5*K2;
                        K4 = dtout*ax + K3;
                        vx = vx_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*vy;
                        K2 = dtout*vy + 0.5*K1;
                        K3 = dtout*vy + 0.5*K2;
                        K4 = dtout*vy + K3;
                        dy = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ay;
                        K2 = dtout*ay + 0.5*K1;
                        K3 = dtout*ay + 0.5*K2;
                        K4 = dtout*ay + K3;
                        vy = vy_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        xpos = xmax;
                        ypos = ymax;
                        vx = -vx;
                        vy = -vy;
                        tempo = tempo - (dt-dtout);
                        
                    }
                }
                else if((testx) < xmin && (testy) < ymin){
                    if(ax == 0.0 && ay == 0.0){ // no acceleration, simple reflection
                        xpos = xmin + Math.abs(testx - xmin); //xpos + dy;
                        ypos = ymin + Math.abs(testy - ymin); //ypos + dy;
                        vx = -vx;
                        vy = -vy;
                    }
                    else{
                        dyout = dy - (testy - ymin );
                        vyout = vy_old + vario*(Ey - vx_old * Bz);
                        dtout = dyout/vyout; // time  boundary
                        
                        ax = charge_value/mass_value*(Ex + (vy_old+vy)*0.5 * Bz);
                        ay = charge_value/mass_value*(Ey - (vx_old+vx)*0.5 * Bz);

                        K1 = dtout*vx;
                        K2 = dtout*vx + 0.5*K1;
                        K3 = dtout*vx + 0.5*K2;
                        K4 = dtout*vx + K3;
                        dx = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ax;
                        K2 = dtout*ax + 0.5*K1;
                        K3 = dtout*ax + 0.5*K2;
                        K4 = dtout*ax + K3;
                        vx = vx_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*vy;
                        K2 = dtout*vy + 0.5*K1;
                        K3 = dtout*vy + 0.5*K2;
                        K4 = dtout*vy + K3;
                        dy = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ay;
                        K2 = dtout*ay + 0.5*K1;
                        K3 = dtout*ay + 0.5*K2;
                        K4 = dtout*ay + K3;
                        vy = vy_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        xpos = xmin;
                        ypos = ymin;
                        tempo = tempo - (dt-dtout);
                        vy = -vy;
                        vx = -vx;
                    }
                }
                else if((testy) < ymin && (testx) > xmax){
                    if(ax == 0.0 && ay == 0.0){ // no acceleration, simple reflection
                        xpos = xmax - (testx - xmax); //xpos - dy;
                        ypos = ymin + Math.abs(testy-ymin); //ypos + dy;
                        vx = -vx;
                        vy = -vy;
                    }
                    else{
                        dyout = dy - (testy - ymin );
                        vyout = vy_old + vario*(Ey - vx_old * Bz);
                        dtout = dyout/vyout; // time  boundary
                        
                        ax = charge_value/mass_value*(Ex + (vy_old+vy)*0.5 * Bz);
                        ay = charge_value/mass_value*(Ey - (vx_old+vx)*0.5 * Bz);

                        K1 = dtout*vx;
                        K2 = dtout*vx + 0.5*K1;
                        K3 = dtout*vx + 0.5*K2;
                        K4 = dtout*vx + K3;
                        dx = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ax;
                        K2 = dtout*ax + 0.5*K1;
                        K3 = dtout*ax + 0.5*K2;
                        K4 = dtout*ax + K3;
                        vx = vx_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*vy;
                        K2 = dtout*vy + 0.5*K1;
                        K3 = dtout*vy + 0.5*K2;
                        K4 = dtout*vy + K3;
                        dy = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ay;
                        K2 = dtout*ay + 0.5*K1;
                        K3 = dtout*ay + 0.5*K2;
                        K4 = dtout*ay + K3;
                        vy = vy_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        xpos = xmax;
                        ypos = ymin;
                        tempo = tempo - (dt-dtout);
                        vy = -vy;
                        vx = -vx;
                        
                    }
                }
                else if((testy) > ymax && (testx) < xmin){
                    if(ax == 0.0 && ay == 0.0){ // no acceleration, simple reflection
                        xpos = xmin + Math.abs(testx - xmin); //xpos + dy;
                        ypos = ymax - (testy-ymax); //ypos - dy;
                        vx = -vx;
                        vy = - vy;
                    }
                    else{
                        dyout = dy - (testy - ymax);
                        vyout = vy_old + vario*(Ey - vx_old * Bz);
                        dtout = dyout/vyout; // time  boundary
                        
                        ax = charge_value/mass_value*(Ex + (vy_old+vy)*0.5 * Bz);
                        ay = charge_value/mass_value*(Ey - (vx_old+vx)*0.5 * Bz);

                        K1 = dtout*vx;
                        K2 = dtout*vx + 0.5*K1;
                        K3 = dtout*vx + 0.5*K2;
                        K4 = dtout*vx + K3;
                        dx = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ax;
                        K2 = dtout*ax + 0.5*K1;
                        K3 = dtout*ax + 0.5*K2;
                        K4 = dtout*ax + K3;
                        vx = vx_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*vy;
                        K2 = dtout*vy + 0.5*K1;
                        K3 = dtout*vy + 0.5*K2;
                        K4 = dtout*vy + K3;
                        dy = 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        K1 = dtout*ay;
                        K2 = dtout*ay + 0.5*K1;
                        K3 = dtout*ay + 0.5*K2;
                        K4 = dtout*ay + K3;
                        vy = vy_old + 1.0/6.0 * (K1+K4) + 1.0/3.0 * (K2 + K3);

                        xpos = xmin;
                        ypos = ymax;
                        vx = -vx;
                        vy = -vy;
                        tempo = tempo - (dt-dtout);
                        
                    }
                }
            }
        }
    }
}