/*$Id: maketernary.c 892 2006-03-10 15:24:59Z rcatwood $*/
/****************************************************************/
/*   Copyright (c) 1998 - 2004 Dept. of Materials, ICSTM        */
/*   All Rights Reserved                                        */
/*   THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF ICSTM       */
/*   The copyright notice above does not evidence any           */
/*   actual or intended publication of such source code,        */
/*   and is an unpublished work by Dept. of Materials, ICSTM.   */
/*   This material contains CONFIDENTIAL INFORMATION that       */
/*   is the property of Imperial College. Any use,              */
/*   duplication or disclosure not specifically authorized      */
/*   by Imperial College is strictly prohibited.                */
/****************************************************************/
/* This code is part of the umats routines developed at in the  */
/* Materials Processing Group, Dept. of Materials, ICSTM.       */
/*      email p.d.lee or r.atwood @ic.ac.uk for details         */
/****************************************************************/
/*This file was created at Thu Jul 15 17:09:31 BST 2004 by rcatwood on hive.beowulf.cluster */
#define MAKETERNARY_C_REV "maketernary.c $Id: maketernary.c 892 2006-03-10 15:24:59Z rcatwood $"
/* Include the common header for this project */
#include <stdio.h>
#include <math.h>
#include "../machine.h"
typedef double FLOAT;

      int n_isteps = 1000;
      FLOAT Ca_istep ;
      FLOAT Cb_istep ;
      FLOAT ia,ib,Ta,Tb;
      FLOAT Cpa,Cpb;
      FLOAT dT;
      int ctr=0,idx;

      int Cas;
      int Cbs;
   FLOAT Tpure=650,m_a=-12,m_b=-6,k_a=0.1,k_b=0.15;
   FLOAT Tliq,Tsol;
   FLOAT Te,Tmin=400,Tmax=660,Tstep=1;
   FLOAT CaMin=0,CaMax=20,CaStep=.1;
   FLOAT CbMin=0,CbMax=20,CbStep=.1;
   FLOAT Ceut_a=12,Ceut_b=6;
   FLOAT Teut_a,Teut_b;
   FLOAT maxeut;
   int nsteps,ntsteps,ncasteps,ncbsteps;
   int slice;
   int ndata;

void draw_region(unsigned char * data,unsigned char  reg,FLOAT t_S_a,FLOAT t_S_b,FLOAT t_C_a,FLOAT t_C_b,FLOAT toff){
      int i,j;
      int Tidx;
      FLOAT Cia,Cib;
      FLOAT Sia,Sib;
      FLOAT Castep,Cbstep,Sastep,Sbstep;
      FLOAT Ts,Tc;
      FLOAT dCa,dCb,dTe;
      FLOAT ofact,ofacti,ofactj;
      FLOAT Ls,Lc;
      FLOAT lp,lps,lpc,ratio,Lab,Lsc,Dca,Dcb;

      Castep = t_C_a / n_isteps;
      Cbstep = t_C_b / n_isteps;
      Sastep = t_S_a / n_isteps;
      Sbstep = t_S_b / n_isteps;

      Ls=hypot(t_S_a ,t_S_b); 
      Lc=hypot(t_C_a,t_C_b);

      dT = Teut_a - Teut_b;

      for(j=0;j<n_isteps;j++){
         ofact = ( 1 - (FLOAT)(abs(j - n_isteps/2))/(n_isteps/2));
         Sia = t_S_a - j * Sastep; 
         Cia = t_C_a - j * Castep; 
         Sib = j * Sbstep; 
         Cib = j * Cbstep; 
         lps = hypot( (j* Sastep),(j*Sbstep));
         lpc = hypot( (j* Castep),(j*Cbstep));
         Ts = Teut_a - (lps/Ls) * dT;
         Tc = Teut_a - (lpc/Lc) * dT + ofact*toff; 
         for(i=0;i<n_isteps;i++){
            
            Lsc = hypot( (Sia-Cia),(Sib-Cib));

            Dca = i * ((Cia-Sia)/n_isteps);
            Cpa = Cia - Dca;
            Dcb = i*((Cib-Sib)/n_isteps);
            Cpb = Cib -Dcb;
            lp = hypot(Dca,Dcb);

            ratio = lp/Lsc;

            Te = Tc - ratio * (Tc - Ts);
            Tidx = (int)(floor( (Te-Tmin) / Tstep));
            //printf("%i Lab,ratio,lp,Cpa,Cpb,Te: %g,%g,%g,%g,%g,%g\n",ctr++,Lab,ratio,lp,Cpa,Cpb,Te);
            Cas = (int)(floor( ((Cpa - CaMin )/ CaStep)));
            Cbs = (int)(floor( ((Cpb - CbMin)/ CbStep)));
            idx= Tidx * slice + Cbs * ncbsteps + Cas;
            if (idx > 0 && idx < ndata && Cas <ncasteps && Cbs < ncbsteps && Tidx < ntsteps ){
               /*
               if (*(data + idx) == 1){
               */
               *(data + idx) = reg;
               /*
               }
               */
            }

         }
      }




}


int main (int argc, char * argv[]){
   unsigned char *data,*datap;
   FLOAT S_a,S_b;
   FLOAT P_a,P_b;

   FLOAT i,j,k;
   FLOAT Ca,Cb;
   FILE * fp;

   ndata = 1;
   S_a = Ceut_a*k_a;
   S_b = Ceut_b*k_b;
   P_a = 86;
   P_b = 70;
   Teut_a = Tpure + m_a * Ceut_a;
   Teut_b = Tpure + m_b * Ceut_b;
   maxeut = MAX(Teut_a,Teut_b);



   ntsteps = (Tmax-Tmin)/Tstep;
   ndata *= ntsteps;

   ncasteps = (CaMax - CaMin)/CaStep;
   ndata *= ncasteps;

   ncbsteps = (CbMax - CbMin)/CbStep;
   ndata *= ncbsteps;
   slice = ncasteps*ncbsteps;
   
   data = (unsigned char *)calloc(ndata,sizeof(unsigned char *));
   datap = data;

   printf("A,B,T: %i %i %i\n",ncasteps,ncbsteps,ntsteps);

   Te=Tmin;
   for (k=0;k<ntsteps;k++){
      Cb=CbMin;
      for(j=0;j<ncbsteps;j++){
         Ca = CaMin;
         for(i=0;i<ncasteps;i++){
         Tliq = Tpure + m_a * Ca + m_b * Cb;
         Tsol = Tpure + (m_a/k_a) * Ca + (m_b/k_b) * Cb;

         if (Te > Tliq) *datap = 0;
         else if (Te> Tsol) *datap = 1;
         else *datap=2;

         datap++;
         Ca += CaStep;
         }
      Cb += CbStep;
      }
   Te += Tstep;
   }

      draw_region(data,5,S_a,S_b,Ceut_a,Ceut_b,30);
      draw_region(data,6,S_a,S_b,P_a,P_b,0);
      draw_region(data,12,P_a,P_b,Ceut_a,Ceut_b,30);

      #ifdef JUNK
      draw_region(data,5,S_a,S_b,Ceut_a,Ceut_b,0);
      draw_region(data,7,Ceut_a,Ceut_b,P_a,P_b,100);
      draw_region(data,9,Ceut_a,Ceut_b,P_a,P_b,0);
      draw_region(data,6,S_a,S_b,P_a,P_b,0);
      Cb=S_b;
      Ca=S_a;
      Ca_istep = (Ceut_a - S_a) / n_isteps;
      Cb_istep = (Ceut_b - S_b) / n_isteps;

      dT = Teut_a - Teut_b;

      for(j=0;j<n_isteps;j++){

         Lab = sqrt(Ca*Ca+Cb*Cb);
         for(i=0;i<n_isteps;i++){
            Dca = i * (Ca/n_isteps);
            Cpa = Ca - Dca;
            Cpb = i*(Cb/n_isteps);
            lp = sqrt(Dca * Dca + Cpb * Cpb);

            ratio = lp/Lab;

            Te = Teut_a - ratio * (Teut_a - Teut_b);
            Ts = (int)(floor( (Te-Tmin) / Tstep));
            //printf("%i Lab,ratio,lp,Cpa,Cpb,Te: %g,%g,%g,%g,%g,%g\n",ctr++,Lab,ratio,lp,Cpa,Cpb,Te);
            Cas = (int)(floor( ((Cpa - CaMin )/ CaStep)));
            Cbs = (int)(floor( ((Cpb - CbMin)/ CbStep)));
            idx= Ts * slice + Cbs * ncbsteps + Cas;
            if (idx > 0 && idx < ndata && Cas <ncasteps && Cbs < ncbsteps && Ts < ntsteps ){
               if (*(data + idx) == 1){
               *(data + idx) = 5+j;
               }
            }

         }
      Ca += Ca_istep;
      Cb += Cb_istep;
      }
      #endif

   fp = fopen("pd.raw","w");
   fwrite(data,sizeof(unsigned char),ndata,fp);
   fclose (fp);
   printf("A,B,T: %i %i %i\n",ncasteps,ncbsteps,ntsteps);
}


