import java.io.*;
import mpi.*;

//30 Apr 2002
//Au-Ni binary phase diagram
//Set Ni A component, Au B component, Solid A phase, Liquid B phase


public class PhaseMPI
{
	public static void main(String[] args) throws MPIException
	{
		PrintStream pfout = null;
		PrintStream iout = null;
		int guest_processor=1, host_processor=0, myrank, worldsize, con_res = 0, h=0;
		double temp = 0, T_STEP, T_MAX, T_START;
		MPI.Init(args);
		myrank = MPI.COMM_WORLD.Rank();
		worldsize  = MPI.COMM_WORLD.Size();
		double timer_array[] = new double [2];
		System.out.println("\n"+myrank+": Started");
		MPI.COMM_WORLD.Barrier();
		timer_array[0] = MPI.Wtime();

		if ( myrank == host_processor )
//I am the host processor
		{ 
			System.out.println("\n\n MPI Started \n");
			try 
			{
				FileOutputStream fout = new FileOutputStream("AuNi1.dat");
				pfout = new PrintStream(fout);
				FileOutputStream infoOut = new FileOutputStream("MPIinfo10.dat");
				iout = new PrintStream(infoOut);
			} 
				catch(IOException ioe) 
			{
				System.err.println("Error in FileOutput...."+ioe);
			}
		}
		//these define the start and stop of the temperature loop in Kelvin
		T_START = 300;
		T_MAX = 1600;
		T_STEP = 1;
		//con_res defines the steps in consintration, the step size across 
		//the concentration axis is 1/con_res
		con_res = 1000;
		for (temp = T_START; temp < T_MAX; temp += T_STEP)
		{
			if (myrank == guest_processor) //Other computers
			{
				CALPHAD auni=new CALPHAD();
				TestG test=new TestG();
				test.setOmega(21689,0);
				test.setT(temp);
				double G0Aa=auni.getG0Aa(temp);
				test.setG0Aa(G0Aa);
				double G0Ab=auni.getG0Ab(temp);
				test.setG0Ab(G0Ab);
				double G0Ba=auni.getG0Ba(temp);
				test.setG0Ba(G0Ba);
				double G0Bb=auni.getG0Bb(temp);
				test.setG0Bb(G0Bb);
				double results[]=test.runTest(con_res);
				MPI.COMM_WORLD.Send(results,0,2*con_res,MPI.DOUBLE,host_processor,10);
			}
			if (myrank == host_processor) //Host Computer
			{
				double receive[] = new double [2*con_res];
				if(temp%10==0){System.out.println(temp);}
				double TdegC=temp-273.0;
				MPI.COMM_WORLD.Recv(receive, 0 ,2*con_res, MPI.DOUBLE, guest_processor, 10);
				while(receive[h]!=0)
				{	
					pfout.println(receive[h] + "\t" + TdegC);
					h++;
				}
				h=0;
			}
			//If we reach the end of the processors, 
			//then reset the current processor to 1 otherwise advance
			if (guest_processor == (worldsize-1))
			{
				guest_processor=1;
			}
			else
			{
				guest_processor++;
			}
		}
		timer_array[1] = MPI.Wtime();
		MPI.Finalize();
		if (myrank==0)
		{
			iout.println("\n\n\t Total Processors: "+worldsize+" (including host),\n");
			iout.println("\t Processors computing: "+ (worldsize-1));
			iout.println("\n\t MPI Total time: "+(timer_array[1]-timer_array[0])+"\n");
		}
	}
}



class TestG
{
	final double R=8.314;
	double G0Aa;
	double G0Ab;
	double G0Ba;
	double G0Bb;
	double T;
	double NA_beta;
	double Omega_a;
	double Omega_b;
	int res;
	
	public TestG()
	{}
	public void setT(double temperature)
	{
		T=temperature;
	}
	public void setG0Aa(double GAa)
	{
		G0Aa=GAa;
	}
	public void setG0Ab(double GAb)
	{
		G0Ab=GAb;
	}
	public void setG0Ba(double GBa)
	{
		G0Ba=GBa;
	}
	public void setG0Bb(double GBb)
	{
		G0Bb=GBb;
	}
	
	public void setNA_beta(double N)
	{
		NA_beta=N;
	}
	
	public void setOmega(double Oa, double Ob)
	{
		Omega_a=Oa;
		Omega_b=Ob;
	}
	public boolean checkA(double NA_alpha)
	{
		double G_a=G0Aa+R*T*Math.log(NA_alpha)
		+Omega_a*Math.pow((1-NA_alpha),2);
		double G_b=G0Ab+R*T*Math.log(NA_beta)
		+Omega_b*Math.pow((1-NA_beta),2);
		double check=G_a-G_b;
		boolean A;
		if((check>(-10)) && (check<(10)))
		{A=true;}
		else {A=false;}
		return A;
	}
	public boolean checkB(double NA_alpha)
	{
		double G_a=G0Ba+R*T*Math.log(1-NA_alpha)
		+Omega_a*Math.pow((NA_alpha),2);
		double G_b=G0Bb+R*T*Math.log((1-NA_beta))
		+Omega_b*Math.pow((NA_beta),2);
		double check=G_a-G_b;
		boolean B;
		if((check>(-10)) && (check<(10)))
		{B=true;}
		else {B=false;}
		return B;
	}
	public double[] runTest(int resolution)
	{
		res=resolution;
		int k=0;
		double internal[] = new double [2*res];
		for (int b=1;b<res;b++)
		{
			double NAb=(double)b/res;
			this.setNA_beta(NAb);
			for (int a=1;a<res;a++)
			{
				double NAa=(double)a/res;
				boolean A=this.checkA(NAa);
				boolean B=this.checkB(NAa);
				if((A==true)&&(B==true))
				{
					
					internal[k] = NAa;
					internal[k+1] = NAb;
					k+=2;
				}
			}
		}
		return internal;
	}
}

/* 
The CRC class is not used in this analsis, although it is an accecptable method
for solving the Gibbs Free Energy.
*/
class CRC
{
	double a,b,c,d,A,B,G;
	
	public CRC(double aa, double bb, double cc, double dd, double AA, 
			double BB)
	{
		a=aa;
		b=bb;
		c=cc;
		d=dd;
		A=AA;
		B=BB;
	}
	
	public double findGcrc(double T)
	{
		G=((-a)*T*Math.log(T))-(0.5*b*0.001*Math.pow(T,2))-
		(0.166666666*c*0.0000001*Math.pow(T,3))-((0.5*d*100000)/T)+
		(B*T+a*T)-(A*1000);
		return G*4.1868;
	}
}

//Must Edit CALPHAD for specific fits.
class CALPHAD
{
	public CALPHAD()
	{}
	public double getG0Aa(double T)
	{
		double G0NiS=0;
		if ((T>=298)&&(T<=1728))
		{
			G0NiS=this.findG_Ni_s_298_1728(T);
		}
		else if(T>1728)
		{
			G0NiS=this.findG_Ni_s_1728_3000(T);
		}
		return G0NiS;
	}
	public double getG0Ba(double T)
	{
		double G0AuS=0;
		if((T>=298)&&(T<=933))
		{
			G0AuS=this.findG_Au_s_298_933(T);
		}
		else if((T>933)&&(T<=1337))
		{
			G0AuS=this.findG_Au_s_933_1337(T);
		}
		else if((T>1337)&&(T<=1735))
		{
			G0AuS=this.findG_Au_s_1337_1735(T);
		}
		else if(T>1735)
		{
			G0AuS=this.findG_Au_s_1735_3200(T);
		}
		return G0AuS;
	}
	public double getG0Ab(double T)
	{
		double G0NiL=0;
		if ((T>=298)&&(T<=1728))
		{
			G0NiL=this.findG_Ni_l_298_1728(T);
		}
		else if(T>1728)
		{
			G0NiL=this.findG_Ni_l_1728_3000(T);
		}
		return G0NiL;
	}
	public double getG0Bb(double T)
	{
		double G0AuL=0;		
		if((T>=298)&&(T<=933))
		{
			G0AuL=this.findG_Au_l_298_933(T);
		}
		else if((T>933)&&(T<=1337))
		{
			G0AuL=this.findG_Au_l_933_1337(T);
		}
		else if((T>1337)&&(T<=1735))
		{
			G0AuL=this.findG_Au_l_1337_1735(T);
		}
		else if(T>1735)
		{
			G0AuL=this.findG_Au_l_1735_3200(T);
		}
		return G0AuL;
	}

	public double findG_Au_l_298_933(double T)
	{
		double G=5613.147
		+97.446385*T
		-22.75455*T*Math.log(T)
		-3.85924*Math.pow(10,-3)*Math.pow(T,2)
		+0.379625*Math.pow(10,-6)*Math.pow(T,3)
		-25097*Math.pow(T,-1);
		return G;
	}
	public double findG_Au_l_933_1337(double T)
	{
		double G=-81023.261
		+1012.217324*T
		-155.6947*T*Math.log(T)
		+87.56015*Math.pow(10,-3)*Math.pow(T,2)
		-11.518713*Math.pow(10,-6)*Math.pow(T,3)
		+10637210*Math.pow(T,-1);
		return G;
	}
	public double findG_Au_l_1337_1735(double T)
	{
		double G=326614.987
		-2025.7579*T
		+263.2523*T*Math.log(T)
		-118.21685*Math.pow(10,-3)*Math.pow(T,2)
		+8.923845*Math.pow(10,-6)*Math.pow(T,3)
		-67999850*Math.pow(T,-1);
		return G;
	}
	public double findG_Au_l_1735_3200(double T)
	{
		double G=413.343
		+155.893158*T
		-30.9616*T*Math.log(T);
		return G;
	}
	public double findG_Au_s_298_933(double T)
	{
		double G=-6938.853
		+106.830495*T
		-22.75455*T*Math.log(T)
		-3.85924*Math.pow(10,-3)*Math.pow(T,2)
		+0.379625*Math.pow(10,-6)*Math.pow(T,3)
		-25097*Math.pow(T,-1);
		return G;
	}
	public double findG_Au_s_933_1337(double T)
	{
		double G=-93575.261
		+1021.601434*T
		-155.6947*T*Math.log(T)
		+87.56015*Math.pow(10,-3)*Math.pow(T,2)
		-11.518713*Math.pow(10,-6)*Math.pow(T,3)
		+10637210*Math.pow(T,-1);
		return G;
	}
	public double findG_Au_s_1337_1735(double T)
	{
		double G=314062.987
		-2016.37379*T
		+263.2523*T*Math.log(T)
		-118.21685*Math.pow(10,-3)*Math.pow(T,2)
		+8.923845*Math.pow(10,-6)*Math.pow(T,3)
		-67999850*Math.pow(T,-1);
		return G;
	}
	public double findG_Au_s_1735_3200(double T)
	{
		double G=-12138.657
		+165.277268*T
		-30.09616*T*Math.log(T);
		return G;
	}
	public double findG_Ni_l_298_1728(double T)
	{
		double G=11235.527
		+108.457*T
		-22.096*T*Math.log(T)
		-4.8407*Math.pow(10,-3)*Math.pow(T,2)
		-382.318*Math.pow(10,-23)*Math.pow(T,7);
		return G;
	}
	public double findG_Ni_l_1728_3000(double T)
	{
		double G=-9549.775
		+268.598*T
		-43.1*T*Math.log(T);
		return G;
	}
	public double findG_Ni_s_298_1728(double T)
	{
		double G=-5179.159
		+117.854*T
		-22.096*T*Math.log(T)
		-4.8407*Math.pow(10,-3)*Math.pow(T,2);
		return G;
	}
	public double findG_Ni_s_1728_3000(double T)
	{
		double G=-27840.655
		+279.135*T
		-43.1*T*Math.log(T)
		+1127.54*Math.pow(10,28)*Math.pow(T,-9);
		return G;
	}
}