
/****************************************************************/
/*      Copyright (c) 1993 Peter D Lee                          */
/*      Copyright (c) 1998 Dept. of Materials, ICSTM            */
/*      All Rights Reserved                                     */
/*      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.*/
/*      continuing D Phil work from University of Oxford        */
/****************************************************************/
/* 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 @imperial.ac.uk for details   */
/****************************************************************/

/********************************************************************************/
/*  This version is distributed under a BSD style public license, as follows:   */
/*                                                                              */
/*  Copyright (c) 2007, Dept. of Materials, Imperial College London             */
/*  All rights reserved.                                                        */
/*  Redistribution and use in source and binary forms, with or without          */
/*  modification, are permitted provided that the following conditions          */
/*  are met:                                                                    */
/*                                                                              */
/*  * Redistributions of source code must retain the above copyright            */
/*  notice, this list of conditions and the following disclaimer.               */
/*                                                                              */
/*  * Redistributions in binary form must reproduce the above                   */
/*  copyright notice, this list of conditions and the following                 */
/*  disclaimer in the documentation and/or other materials provided             */
/*  with the distribution.                                                      */
/*                                                                              */
/*  * Neither the name of the Dept. of Materials, Imperial College London, nor  */
/*  the names of its contributors may be used to endorse or promote products    */
/*  derived from this software without specific prior written permission.       */
/*                                                                              */
/*  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS         */
/*  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT           */
/*  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR       */
/*  A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT        */
/*  OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,       */
/*  SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED    */
/*  TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR      */
/*  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF      */
/*  LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING        */
/*  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS          */
/*  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                */
/********************************************************************************/
/*END of LICENSE NOTICE*/

/****************************************************************/
/* READPHASE.C:   (Part of CA code)                             */
/* Subroutine to read the phase diagram from a file.            */
/****************************************************************/
/****************************************************************/
/* Written by Peter D. Lee, Robert C. Atwood and A. Chirazi     */
/*                                            Imperial College */
/* Feb 20, 2001                                                 */
/****************************************************************/

/****************************************************************/
/* Versions maintained with RCS                                 */
/* Version 1.0: Aug. 13, pdl                                    */
/****************************************************************/
/*RCS id $Id: readphase.c 1339 2008-07-23 13:58:29Z  $*/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "machine.h"
#include "read_ctrl.h"
#include "readmat.h"
#include "blocks.h"
#include "multi_diff_props.h"

/*******sort an array of n elements using the SHELL sort method******/
void shell_ca (int n, CA_FLOAT * a)
{
  int i, j, inc;
  CA_FLOAT v;

  inc = 1;

  do {
    inc *= 3;
    inc++;
  } while (inc <= n);
  do {
    inc /= 3;
    for (i = inc + 1; i <= n; i++) {
      v = a[i];
      j = i;
      while (a[j - inc] > v) {
        a[j] = a[j - inc];
        j -= inc;
        if (j <= inc)
          break;
      }
      a[j] = v;
    }
  } while (inc > 1);
}

/*****end of shell_ca routine for sort******/

/*****function to generate tabulated values for bi-linear inter***/
void generate_matrix (Ctrl_str * cp, BB_struct * bp)
{
  int i, j, k, l, h;
  int numtietri, counter;
  int ele_num, ele_1, test;
  MultiS_struct *ms;

  ms = &(bp->MultiSvals);
  ele_num = cp->NUM_COMP;
  ele_1 = ele_num - 1;
  numtietri = ms->numtietri;

 /***********************allocate arrays*****************/

  ms->xxa = (CA_FLOAT **) calloc (ele_1 + 3, sizeof (CA_FLOAT *));

  for (i = 0; i <= ele_1 + 2; i++) {
    ms->xxa[i] = (CA_FLOAT *) calloc (numtietri + 1, sizeof (CA_FLOAT));
  }

 /*************************************************************/

 /**transform the phase diagram into three separate matrices**/

  for (j = 0; j < ele_1 + 2; j++) {
    for (i = 0, counter = 1; i < numtietri; i++) {
      if (i == 0) {
        ms->xxa[j + 1][counter] = ms->phaval[0][j];
      } else {
        for (k = 1, test = 0; k <= counter; k++) {
          if (ms->phaval[i][j] != ms->xxa[j + 1][k]) {
            test++;
          }
        }
        if (test == counter) {
          counter++;
          ms->xxa[j + 1][counter] = ms->phaval[i][j];
        }
      }                         /* end of main if */
    }                           /* end of i loop */
    ms->counter[j + 1] = counter;
  }                             /* end of j loop */

 /*****************begin the matrice generation*************/

 /****************do sorts before matrice generation*******/

  for (j = 1; j <= ele_1 + 2; j++) {
    shell_ca (ms->counter[j], &(ms->xxa[j][0]));
  }
 /****************sort finished****************************/

 /****************allocate arrays for function values*******/

  ms->ya1 = (CA_FLOAT **) calloc (ms->counter[3] + 1, sizeof (CA_FLOAT *));
  ms->ya2 = (CA_FLOAT **) calloc (ms->counter[1] + 1, sizeof (CA_FLOAT *));
  ms->ya3 = (CA_FLOAT **) calloc (ms->counter[1] + 1, sizeof (CA_FLOAT *));

  for (i = 0; i <= ms->counter[3]; i++) {
    ms->ya1[i] = (CA_FLOAT *) calloc (ms->counter[4] + 1, sizeof (CA_FLOAT));
  }
  for (i = 0; i <= ms->counter[1]; i++) {
    ms->ya2[i] = (CA_FLOAT *) calloc (ms->counter[4] + 1, sizeof (CA_FLOAT));
  }
  for (i = 0; i <= ms->counter[1]; i++) {
    ms->ya3[i] = (CA_FLOAT *) calloc (ms->counter[3] + 1, sizeof (CA_FLOAT));
  }

 /*************************************************************/

  for (i = 0; i < numtietri; i++) {
    for (j = 1; j <= ms->counter[3]; j++) {
      for (k = 1; k <= ms->counter[4]; k++) {
        if (ms->xxa[3][j] == ms->phaval[i][2] && ms->xxa[4][k] == ms->phaval[i][3]) {
          ms->ya1[j][k] = ms->phaval[i][0];
        }
      }
    }
  }

  for (i = 0; i < numtietri; i++) {
    for (j = 1; j <= ms->counter[1]; j++) {
      for (k = 1; k <= ms->counter[4]; k++) {
        if (ms->xxa[1][j] == ms->phaval[i][0] && ms->xxa[4][k] == ms->phaval[i][3]) {
          ms->ya2[j][k] = ms->phaval[i][2];
        }
      }
    }
  }

  for (i = 0; i < numtietri; i++) {
    for (j = 1; j <= ms->counter[1]; j++) {
      for (k = 1; k <= ms->counter[3]; k++) {
        if (ms->xxa[1][j] == ms->phaval[i][0] && ms->xxa[3][k] == ms->phaval[i][2]) {
          ms->ya3[j][k] = ms->phaval[i][3];
        }
      }
    }
  }

 /**************************generation finished*************************/

}

/******************************end of generate_matrix********************/

int readphase (Ctrl_str * cp, BB_struct * bp)
{
  char line[MAX_STRING_LEN];
  char *token;
  int i, j, itmp, n, npaded;
  int nx, ny, nz, nc;
  int rflag = 0;
  int index = 0;
  int error = FALSE;
  Mat_str *mp;
  Nuc_str *np;
  P_str *pp;
  MultiS_struct *ms;
  int ele_num, ele_1;
  CA_FLOAT *Cinit_multi_ptr;
  CA_FLOAT **Clim_multi_ptr;
  FILE *fileout;

  ms = &(bp->MultiSvals);       /* local variable for the multi diff values */
  pp = &(bp->pprops);
  mp = &(bp->mprops);
  np = &(bp->nprops);

  nx = bp->nc[0];
  ny = bp->nc[1];
  nz = bp->nc[2];
  nc = bp->ncsb;
  npaded = (nx + 2) * (ny + 2) * (nz + 2);
  bp->ntsb = bp->nsb[0] * bp->nsb[1] * bp->nsb[2];

  /*from read_ctrl.h */
  ele_num = D_NUM_COMP;
   /*******/
  ele_num = cp->NUM_COMP;       /*number of components in the alloy */
  ele_1 = ele_num - 1;

  fileout = fopen ("phase_out", "a+");

/*******************read liquidus temperature and solute concentration****/
/*******************from the input file provided by thermocalc************/

   /*********************************************************/
  /* Open the file                                         */
   /*********************************************************/
  if ((cp->fd_phadia = fopen (cp->fn_phadia, "r")) == NULL) {
    fprintf (stderr, "Error: can't open input file [%s]\n", cp->fn_phadia);
    exit (0);
  }

  while (fgets (line, MAX_STRING_LEN, cp->fd_phadia) != NULL) {
    /* ignore comment and blank lines */
    if (line[0] == '%' || line[0] == '#' || line[0] == '\n' || (token = strtok (line, " ,;\t")) == NULL) {
      continue;

      /*********************************************************/
      /* values for liquidus and solute concentration in liquid */
      /*********************************************************/
      /* number of tie-triangles produced by thermocalc */
    } else if (strcasecmp (token, "NumTieTri") == 0) {
      if ((token = strtok (NULL, " ,;\t")) != NULL) {
        ms->numtietri = atoi (token);
      } else {
        ms->numtietri = D_NUMTIETRI;
      }
    } else {
            /*****assign the first token for each loop****/
      ms->phaval[index][0] = atof (token);
      for (i = 1; i <= 3; i++) {
        if ((token = strtok (NULL, " ,;\t")) != NULL) {
          ms->phaval[index][i] = atof (token);
        }
      }
      index++;
    }
  }                             /* end of while routine */

/******************call generate_matrix routine********************/

  generate_matrix (cp, bp);

/*******************print test for the phase diagram routine******/

  for (i = 0; i < ms->numtietri; i++) {
    fprintf (fileout, "%d \t %f \t %f \t %f \t %f \n", i, ms->phaval[i][0], ms->phaval[i][1], ms->phaval[i][2], ms->phaval[i][3]);
  }
  for (i = 1; i <= 4; i++) {
    fprintf (fileout, "row number %d with %d elements is printed below \n \n", i, ms->counter[i]);
    for (j = 1; j <= ms->counter[i]; j++) {
      fprintf (fileout, "%d \t %f \n", j, ms->xxa[i][j]);
    }
  }
  for (i = 1; i <= ms->counter[3]; i++) {
    for (j = 1; j <= ms->counter[4]; j++) {
      fprintf (fileout, "liquidus function, index %d and %d = %f for %f and %f  \n", i, j, ms->ya1[i][j], ms->xxa[3][i],
               ms->xxa[4][j]);
    }
  }
  fclose (fileout);
  return (1);

}

/***end of readphase subroutine***/
/* Little subroutine to get rcs id into the object code */
/* so you can use ident on the compiled program  */
/* also you can call this to print out or include the rcs id in a file*/
char const *rcs_id_readphase_c ()
{
  static char const rcsid[] = "$Id: readphase.c 1339 2008-07-23 13:58:29Z  $";

  return (rcsid);
}

/* end of rcs_id_readphase_c subroutine */
