/****************************************************************************
 ***                                                                      ***
 ***      S T A R   P R O B E        F I T T E R                          ***
 ***                                                                      ***
 ****************************************************************************/
/*****
 ***  This module evaluates how far the inward & outward integrations
 ***  are from each other, and adjust the boundary values of (p,t,l,r)
 ***  accordingly.
 *****/

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

extern FILE *pf;

double delta     [4],
       old_delta [4],
       a         [4] [8];

/*****
 ***  D i s c o n t i n u i t y
 ***
 ***  Returns (true) if any of the endpoints are more than <epsilon> apart.
 *****/

int discontinuity()
{
 iam(70);
 int i,count,val;
 double epsilon;
 block(in,"discontinuity",0.0);
 
 epsilon = 1.0e-2;

 count = 0;
 for (i=0;i<4;i++){
   old_delta[i] = delta[i];
   delta[i] = results[i][0] - results[i][1];
   if ((abs(delta[i]/results[i][0])) > epsilon)
     count++;
   }

 block(mid,"discontinuity count",(double)count);
 val =  (count > 0) ? 1 : 0;

 block(out,"discontinuity",(double)val);
 return(val);
}

/*****
 ***  X f e r   R e s u l t s
 ***
 ***  Transfers results array into 'a[]' 
 *****/

void xfer_results()
{
 iam(71);
 register int i;
 double e1,e2,e3,e4;
 block(in,"xfer results",0.0);

 e1 = corepress    * (adj_p - 1.0);
 e2 = coretemp     * (adj_t - 1.0);
 e3 = lumratio*SL  * (adj_l - 1.0);
 e4 = radratio*SR  * (adj_r - 1.0);

 for (i=0;i<4;i++){
   a[i][0] = (results[i][2] - results[i][0])/e1;
   a[i][1] = (results[i][3] - results[i][0])/e2;
   a[i][2] = (results[i][4] - results[i][1])/e3;
   a[i][3] = (results[i][5] - results[i][1])/e4;
   }

 block(out,"xfer results",0.0);
}

/*****
 ***  I n v e r t
 ***
 ***  Inverts the 4x4 matrix a[] into a'[]
 *****/

void dump_a()
{
int i,j;
 fprintf(pf,"\naaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa");
 for (i=0;i<4;i++){
   fprintf(pf,"\na");
   for (j=0;j<8;j++)
     fprintf(pf," %g",a[i][j]);
   }
 fprintf(pf,"\naaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa\n");
}

void invert()
{
 iam(72);
 register int i,j,im,ip,ist,n,n1;
 double al;
 block(in,"invert",0.0);

for (i=0;i<4;i++)
   for (j=4;j<8;j++)
     a[i][j] = (i == (j-4)) ? 1.0 : 0.0;

 if (trace_me) dump_a();

 n1 = 8; n = 4;

 for (ip=0;ip<n;ip++){
   im = ip;
   ist = ip + 1;
   if (ip != n)
     for (i=ist;i<n;i++)
       if (abs(a[im][ip]) < abs(a[i][ip]))
         im = i;
   if (a[im][ip] == 0.0)
     assert_error(proc_num,(double)im,(double)ip,(double)ist,0.0,0.0,0.0,0.0);
   if (im != ip)
     for (j=ip;j<n1;j++){
       al=a[ip][j]; a[ip][j]=a[im][j]; a[im][j]=al;}
   al = a[ip][ip];
   a[ip][ip] = 1.0;
   for (j=ist;j<n1;j++)
     a[ip][j] /= al;
   for (i=0;i<n;i++)
     if (i!=ip){
       al = a[i][ip];
       for (j=ip;j<n1;j++)
         a[i][j] -= al*a[ip][j];
       }
   }
 if (trace_me) dump_a();
 block(out,"invert",0.0);
}
 

/*****
 ***  N e w   B o u n d a r y   C o n d i t i o n s
 ***
 ***  Ref: Sears & Brownlee (Aller) pp. 597-601
 *****/

void new_boundary_conditions()
{
 iam(73);
 register int i,j,k;
 double new_val[4],dp,dt,dl,dr;
 block(in,"new bdry cond",0.0);

 xfer_results();

 invert();

 for (i=0;i<4;i++){
  new_val[i] = 0.0;
  for (j=4,k=0;j<8;j++,k++)
    new_val[i] += (a[i][j] * delta[k]);
  }

 dp = new_val[0]/corepress;
 dt = new_val[1]/coretemp;
 dl = new_val[2]/(SL*lumratio);
 dr = new_val[3]/(SR*radratio);

 block(mid,"new bdry cond +p",new_val[0]);
 block(mid,"new bdry cond +t",new_val[1]);
 block(mid,"new bdry cond +l",new_val[2]);
 block(mid,"new bdry cond +r",new_val[3]);

 printf("... new bdry cond: %g %g %g %g\n",dp,dt,dl,dr);
 fprintf(pf,"... new bdry cond: %g %g %g %g\n",dp,dt,dl,dr);

 corepress += new_val[0];
 coretemp  += new_val[1];
 lumratio  = ((SL*lumratio)+new_val[2])/SL;
 radratio  = ((SR*radratio)+new_val[3])/SR;

 block(out,"new bdry cond",0.0);
}


/*****
 ***  C l o s e r   F i t
 ***
 ***  Returns (true) if all (p,t,l,r) deltas were reduced from the previous
 ***  base integrations (0 & 1).
 *****/

int closer_fit(ixstop)
int ixstop;
{
 iam(74);
 int val;
 double dP,dT,dL,dR;
 block(in,"closer fit",0.0);

 dP = abs(pressgrad[ixstop]-pressgrad[ixstop+1])/
      ((pressgrad[ixstop]+pressgrad[ixstop+1])*0.5);
 dT = abs(tempgrad[ixstop]-tempgrad[ixstop+1])/
      ((tempgrad[ixstop]+tempgrad[ixstop+1])*0.5);
 dL = abs(lumgrad[ixstop]-lumgrad[ixstop+1])/
      ((lumgrad[ixstop]+lumgrad[ixstop+1])*0.5);
 dR = abs(radiusgrad[ixstop]-radiusgrad[ixstop+1])/
      ((radiusgrad[ixstop]+radiusgrad[ixstop+1])*0.5);

 old_diff = cur_diff;
 cur_diff = (dP+dT+dL+dR)*0.25;

 block(mid,"closer fit dP",dP);
 block(mid,"closer fit dT",dT);
 block(mid,"closer fit dL",dL);
 block(mid,"closer fit dR",dR);

 val = (cur_diff < old_diff) ? 1 : 0;

 block(out,"closer fit",(double)val);
 return(val);
}