/****************************************************************************
 ***                                                                      ***
 ***      S T A R   P R O B E       I N I T I A L I Z E                   ***
 ***                                                                      ***
 ****************************************************************************/
/*****
 ***  This module opens libraries, opens files, gets options, and
 ***  gets initial parameters.
 *****/

#include "stdio.h"
#include "math.h"
#include "limits.h"
#include <dos.h>
#include "spmac.h"
#include "spref.h"

extern FILE *pf;

/*****
 *** I n i t i a l i z e
 *****/

void initialize()
{
iam(1);
int i,j;
FILE *inf;
double d;

 pf = fopen("sprpt.rpt","w");
 inf= fopen("spparm.in","r");

 printf("**************************************************\n");
 printf("***        S  T  A  R  P  R  O  B  E           ***\n");
 printf("***                                            ***\n");
 printf("***        version 1.0    19-nov-86            ***\n");
 printf("**************************************************\n");
 fprintf(pf,"**************************************************\n");
 fprintf(pf,"***        S  T  A  R  P  R  O  B  E           ***\n");
 fprintf(pf,"***                                            ***\n");
 fprintf(pf,"***        version 1.0    19-nov-86            ***\n");
 fprintf(pf,"**************************************************\n");


/* set defaults */

 for (i=0;i<MAX_PROCS;i++) proc_trace[i]=0;
 massratio = 1.0;
 debug = 0;
 max_tries = 2;
 X = 0.70;  Y = 0.27;  Z = 0.03;
 fitptmass = 0.5;
 massgrad[ 0] = 0.001;
 massgrad[ 1] = 0.001;
 massgrad[ 2] = 0.001;
 massgrad[ 3] = 0.001;
 massgrad[ 4] = 0.001;
 massgrad[ 5] = 0.001;
 massgrad[ 6] = 0.001;
 massgrad[ 7] = 0.001;
 massgrad[ 8] = 0.001;
 massgrad[ 9] = 0.001;
 massgrad[10] = 0.01;
 massgrad[11] = 0.01;
 massgrad[12] = 0.01;
 massgrad[13] = 0.01;
 massgrad[14] = 0.01;
 massgrad[15] = 0.01;
 massgrad[16] = 0.01;
 massgrad[17] = 0.01;
 massgrad[18] = 0.01;
 massgrad[19] = 0.05;
 massgrad[20] = 0.05;
 massgrad[21] = 0.10;
 massgrad[22] = 0.10;
 massgrad[23] = 0.10;
 massgrad[24] = 0.10;
 massgrad[25] = 0.10;
 massgrad[26] = 0.10;
 massgrad[27] = 0.05;
 massgrad[28] = 0.05;
 massgrad[29] = 0.01;
 massgrad[30] = 0.01;
 massgrad[31] = 0.01;
 massgrad[32] = 0.01;
 massgrad[33] = 0.01;
 massgrad[34] = 0.01;
 massgrad[35] = 0.01;
 massgrad[36] = 0.01;
 massgrad[37] = 0.01;
 massgrad[38] = 0.01;
/*
 massgrad[38] = 0.001;
 massgrad[39] = 0.001;
 massgrad[40] = 0.001;
 massgrad[41] = 0.001;
 massgrad[42] = 0.001;
 massgrad[43] = 0.001;
 massgrad[44] = 0.001;
 massgrad[45] = 0.001;
 massgrad[46] = 0.001;
 massgrad[47] = 0.001;
*/
 gradsteps = 39;

/* get parmfile data */

 if (fscanf(inf,"%lf ",&massratio)==EOF) goto EndInit;

 if (fscanf(inf,"%d ",&debug)==EOF) goto EndInit;

 if (debug==1){
   i = 99;
   while (i > 0){
     if (fscanf(inf,"%d ",&i)==EOF) goto EndInit;
     if ((i > 0) & (i < MAX_PROCS)) proc_trace[i] = 1;
     }
   }

 if (debug>0) printf("Debug = %d\n",debug);

 block(in,"Initialize",0.0);
 
 if (fscanf(inf,"%d ",&max_tries)==EOF) goto EndInit;
 
 if (fscanf(inf,"%lf ",&X)==EOF) goto EndInit;
 if (fscanf(inf,"%lf ",&Y)==EOF) goto EndInit;
 if (fscanf(inf,"%lf ",&Z)==EOF) goto EndInit;

 if (fscanf(inf,"%lf ",&fitptmass)==EOF) goto EndInit;
 if (fscanf(inf,"%d ",&gradsteps)==EOF) goto EndInit;
 d=0.0;
 for (i=0;i<gradsteps;i++)
   if (fscanf(inf,"%lf ",&massgrad[i])==EOF){
     printf("*** Insufficient mass elements = %d",i);
     exit(1);
     }
   else d+=massgrad[i];


EndInit:
 mass=SM*massratio;
 u=meanmolweight(X,Y,Z);

 printf("Massratio = %g\n",massratio);
 printf("Trials    = %d\n",max_tries);
 printf("Gradsteps = %d\n",gradsteps);
 printf("FitPtMass = %g\n",fitptmass);
 printf("Total mass fraction = %g\n",d);
 printf("X=%g, Y=%g, Z=%g\n",X,Y,Z);
 if (debug > 0) fprintf(pf,"\n Debug trace on. Level=%d.\n",debug);

 for (i=0;i<4;i++)
   for (j=0;j<6;j++)
     results[i][j] = 0.0;

 adj_p = 1.001;
 adj_t = 1.001;
 adj_l = 1.001;
 adj_r = 1.001;
 cur_diff = 0.0;
 old_diff = 0.0;

 fclose(inf);

 block(out,"Initialize",mass);
}

/*****
 *** S e t   B o u n d a r y   C o n d i t i o n s
 ***
 *** Ref: Bowers Eq. 9.9 - 9.23 
 *****/

void set_boundary_conditions()
{
 iam(2);
 block(in,"set_boundary_conditions",massratio);

 effpress = 6.0e9;
 if (massratio == 1.0){ /* special test case */
   radratio = 0.94;
   lumratio = 0.72;
   coredensity = 90.0;
   coretemp = 13.7e6;
   efftemp = 5900.0;
   effpress = 6.0e9;
   }
 else if (massratio<=2.0){/*p-p cycle*/
   radratio=0.312*pow(u,-0.538)*pow(massratio,0.0769);
   lumratio=49.1*pow(u,7.77)*pow(massratio,5.46);
   coretemp=3.05e7 * pow(u,1.54)*pow(massratio,0.923);
   coredensity=86.0*pow(u,-1.615)*pow(massratio,0.769);
   efftemp=2.73e4 * pow(u,2.21)*pow(massratio,1.327);
   }
 else if (massratio<=3.0){/*CNO cycle*/
   radratio=0.451*pow(u,0.395)*pow(massratio,0.697);
   lumratio=43.5*pow(u,7.30)*pow(massratio,5.18);
   coretemp=1.98e7 * pow(u,0.606)*pow(massratio,0.364);
   coredensity=65.8*pow(u,-0.455)*pow(massratio,-0.909);
   efftemp=2.34e4 * pow(u,1.63)*pow(massratio,0.871);
   }
 else{/*electron scattering*/
   radratio=0.454*pow(u,0.588)*pow(massratio,0.765);
   lumratio=112.0*pow(u,4.0)*pow(massratio,3.0);
   coretemp=2.12e7 * pow(u,0.412)*pow(massratio,0.235);
   coredensity=60.3*pow(u,-1.765)*pow(massratio,-1.294);
   efftemp=2.77e4 * pow(u,0.706)*pow(massratio,0.368);
   }
   /* pressure from Clayton eq 2.10 */
 corepress=N0*K*coredensity*coretemp/u;
 B = 0.8;      /*ratio of gas pressure to total pressure */
 poly = 1.25;  /* polytropic index */

 block(out,"set_boundary_conditions",u);
}


/*****
 *** S h u t d o w n
 *****/

void shutdown()
{
iam(99);
int i;
 block(in,"shutdown",0.0);

 report_results();

 for (i=0;i<MAX_PROCS;i++)
   if (proc_count[i] > 0)
    fprintf(pf,"proc %d = %d\n",i,proc_count[i]);

 block(out,"shutdown",0.0);

 fclose(pf);
}