/* string.java
	Program to calculate the propagation of waves on a string
	Written by Justin Elser 3/3/06
	Final draft written 3/18/06
	*/
	
import java.io.*;

public class string
{
	 public static void main(String[] argv)throws IOException, FileNotFoundException
    {
	    int i, j, k;
	    int length=100;	// length of string
	    int N=length;
	    int maxtime=2800;	// number of time steps to calculate
	    double kappa=0.5;	// friction coefficient
	    double rho=0.1;	// density per unit length
	    double T0=50.0;	// Tension (constant)
	    double g=10.0;	// gravity
	    double D=T0/(rho*g);	// constant
	    double c=Math.sqrt(T0/rho);		// velocity of wave
	    double dx=0.01;		// x-change
	    double dt=0.00015;		// time change
	    double c1=dx/dt;		// delta x /delta t
	    double x;			// x coordinate
	    int switchcase=6;	// 1 for single pluck, 2 for double, 3 for sine, 4 for half sine, 5 for cosh, 6 for travel on caternary
	    
	    
	    double wave[][]=new double [length+1][3];
	    double cosh[]=new double [length+1];
	    
	    PrintWriter data = new PrintWriter(new FileOutputStream("wave.dat"), true);	// file for 3d gnuplot
	    
	    PrintWriter init = new PrintWriter(new FileOutputStream("initial.dat"), true);	// file for 2d gnuplot of initial
	    
	    PrintWriter batch = new PrintWriter(new FileOutputStream("batch"), true);// batch file for gnuplot animation
	    
	    batch.println("set term png");		// set output for gnuplot to png, then use convert to convert to animated gif in a shell
	    
	
	switch(switchcase)
	{
	    
	  case 1:
	    for(i=0;i<=length;i++)	// single pluck
	    {
		    x=((double)(i))/((double)(length));
		    cosh[i]=0.0;	// no gravity
		    
		    if (i<(int)(Math.round(0.8*length)))
			    wave[i][0]=1.25*x;
		    else
			    wave[i][0]=5.0-5.0*x;
	    }
	    break;
	    
	  case 2:
	    for(i=0;i<=length;i++)	// double pluck
	    {
		    x=((double)(i))/((double)(length));
		    cosh[i]=0.0;	// no gravity
		    
		    
		    if (x<0.1)
			    wave[i][0]=0.0;
		    else if (x>=0.1 && x<0.2)
			    wave[i][0]=(10.0*x-1.0);
		    else if (x>=0.2 && x<0.3)
			    wave[i][0]=(-10.0*x+3.0);
		    else if (x>=0.3 && x<0.7)
			    wave[i][0]=0.0;
		    else if (x>=0.7 && x<0.8)
			    wave[i][0]=(10.0*x-7.0);
		    else if (x>=0.8 && x<0.9)
			    wave[i][0]=(-10.0*x+9.0);
		    else
			    wave[i][0]=0.0;
		 
	    }
	    break;
	    
	  case 3:
	    for(i=0;i<=length;i++)	// sine function start
	    {
		    x=((double)(i))/((double)(length));
		    wave[i][0]=Math.sin(2*Math.PI*x);
		    cosh[i]=0.0;	// no gravity
	    }
	    break;
	    
	  case 4:
	    for(i=0;i<=length;i++)	// half sine function start
	    {
		    x=((double)(i))/((double)(length));
		    wave[i][0]=Math.sin(Math.PI*x);
		    cosh[i]=0.0;	// no gravity
	    }
	    break;
	  
	  case 5:
	    for(i=0;i<=length;i++)	// cosh start
	    {
		    wave[i][0]=77.154-D*Math.cosh((i-N/2.)/D);
		    cosh[i]=D/2.*Math.cosh((i-N/2.)/D);
	    }
	    break;
	    
	  case 6:			// traveling sine wave for catenary
	    for(i=0;i<=length;i++)
	    {
		    x=((double)(i))/((double)(length));
		    cosh[i]=D/2.*Math.cosh((i-N/2.)/D);
		    
		    if(x<0.1)
			    wave[i][0]=5*Math.sin(10*Math.PI*x);
		    else
			    wave[i][0]=0.0;
	    }
	    break;
	    
	}	// end of switch statement
	
	for(i=0;i<=length;i++)	// print initial config
	    {
		   init.println(wave[i][0]);
	    }
	    
	    
	    for(i=1;i<length;i++)	// set first time step using that dy/dt|t=0=0 (eqn 17.17 in notes)
	    {
		    x=((double)(i))/((double)(length));
		    //wave[i][1]=wave[i][0]+(c*c/(c1*c1))*(wave[i+1][0]+wave[i-1][0]-2*wave[i][0]);	// no gravity
	    	    wave[i][1]=wave[i][0]+dt*dt*(T0/D*(Math.sinh((i-N/2.)/D)/dx*(wave[i][0]-wave[i-1][0])+T0/(dx*dx)*Math.cosh((i-N/2.)/D)*(wave[i+1][0]+wave[i-1][0]-2*wave[i][0])));
		    // above is for with gravity (friction does not change this equation)
	    }
	    
	    
	    for(j=0;j<=maxtime;j++)	// find the next set of coordinates (eqn 17.14 in notes)
	    {
		    
		    
		    
		    for(i=1;i<length;i++)
		    {
			    x=((double)(i))/((double)(length));
/* no grav or fric */	    //wave[i][2]=2*wave[i][1]-wave[i][0]+(c*c/(c1*c1))*(wave[i+1][1]+wave[i-1][1]-2*wave[i][1]); 
/* no fric */		    //wave[i][2]=2*wave[i][1]-wave[i][0]+dt*dt*(T0/D*(Math.sinh((i-N/2.)/D)/dx*(wave[i+1][1]+wave[i-1][1]-2*wave[i][1])+T0/(dx*dx)*Math.cosh((i-N/2.)/D)*(wave[i+1][1]+wave[i-1][1]-2*wave[i][1])));
/* no grav */		    //wave[i][2]=2*wave[i][1]-wave[i][0]+(c*c/(c1*c1))*(wave[i+1][1]+wave[i-1][1]-2*wave[i][1])-2*kappa*dt/rho*(wave[i][1]-wave[i][0]);
/* everthing */		    wave[i][2]=2*wave[i][1]-wave[i][0]+dt*dt*(T0/D*(Math.sinh((i-N/2.)/D)/dx*(wave[i+1][1]+wave[i-1][1]-2*wave[i][1])+T0/(dx*dx)*Math.cosh((i-N/2.)/D)*(wave[i+1][1]+wave[i-1][1]-2*wave[i][1])))-2*kappa*dt/rho*(wave[i][1]-wave[i][0]);
		    }
		    
		    for(i=0;i<=length;i++)	// final becomes initial for next step
		    {
			    wave[i][0]=wave[i][1];
			    wave[i][1]=wave[i][2];
		    }
		    
		    k=j/30;		// only print out every 30 so that graph isn't cluttered
		    if(j%30==0)		
		    {
			    PrintWriter anim = new PrintWriter(new FileOutputStream("anim"+k+".dat"), true);
			    
			    for(i=0;i<=length;i++)
			    {
				    data.println(wave[i][2]+cosh[i]);	// store data in file for 3d plot
				    anim.println(wave[i][2]+cosh[i]);	// store data in files for animation
			    }
			    data.println(" "); 		// insert space for gnuplot
			    if(k<10)
			    {
				    batch.println("set output 'anim00"+k+".png'");
				    batch.println("plot [0:100] [10:50] 'anim"+k+".dat' w l");  // the second [xx:xx] may have to be changed so that plots show up
				    //batch.println("plot 'anim"+k+".dat' w l");
			    }
			    else if(k>10 && k<100)
			    {
				    batch.println("set output 'anim0"+k+".png'");
				    batch.println("plot [0:100] [10:50] 'anim"+k+".dat' w l");
				    //batch.println("plot 'anim"+k+".dat' w l");
			    }
			    else if(k>100)
				    batch.println("set output 'anim"+k+".png'");
				    batch.println("plot [0:100] [10:50] 'anim"+k+".dat' w l");
				    //batch.println("plot 'anim"+k+".dat' w l");
		    }
		    
	    }
    }
}