//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); } } } } } }