/******************************************************************
 * FUNCTION: adiab - Stirling cycle machine simulation            *
 *                   Ideal Adiabatic analysis                     *
 * FILENAME adiab.c                                               *
 * PROGRAMMER I.Urieli    (FORTRAN)                               *
 *            E. Malroy   (C Translation)                         *
 * DATE:     11/30/94  (FORTRAN), 03/12/96 (C Translation)        *
 * LAST MODIFIED: 01/01/97 by I.Urieli                            *
 * GLOBAL VARIABLES: none                                         *
 * SYMBOLIC CONSTANTS: COL,ROWV, ROWD, indices of var matrix      *
 * PROTOTYPE: void adiab(double var[][COL], double dvar[][COL]);  *
 * PRE:   all global variables in define.h have values            *
 * POST:  var, dvar matrices have been filled in                  *
 ******************************************************************/
#include <stdio.h>
#include <math.h>
#include "define.h"
#include "adiabatic.h"
#include "odes.h"

void adiab(double var[][COL],  /* matrix of variables            */
           double dvar[][COL]) /* matrix of derivatives          */
{


   const double epsilon = 1.0; /* allowable error in temperature (K) */
   const int maxit = 20; /* maximum number of iterations to converge */

/*    variables  */

   int i,j, iter, ninc, step, kk;
   double y[ROWV];   /* dependent variable vector */
   double dy[ROWD];    /* dependent derivative vector */
   double tc0;       /* initial cycle value comp temp */
   double te0;       /* initial cycle value exp temp  */
   double terror;    /* cyclic accumulative temp error */
   double theta;     /* angle - independent variable   */
   double dtheta;    /* angle - independent derivative */

   printf("\n...now doing an Ideal Adiabatic simulation...\n");
   fprintf(printfile,"\nIdeal Adiabatic simulation:\n");
   ninc = 360;
   step = ninc/36;
   dtheta = 2.0*PI/ninc;
   y[THE] = th;
   y[TCK] = tk;
   y[TE] = th;
   y[TC] = tk;
   iter = 1;
/*     *** solve **************************************************/
   do {
      printf("\niteration %2d: initial Tc = %.1f(K), Te = %.1f(K)\n",
           iter,y[TC],y[TE]);
      fprintf(printfile,"\niteration %2d: initial Tc = %.1f(K), Te = %.1f(K)\n",
           iter,y[TC],y[TE]);
      tc0 = y[TC];
      te0 = y[TE];
      theta = 0.0;
      y[QK] = 0.0;
      y[QR] = 0.0;
      y[QH] = 0.0;
      y[WC] = 0.0;
      y[WE] = 0.0;
      y[W] = 0.0;
      for(i=0;i<ninc;i++) {
         rk4(dadiab,7,&theta,dtheta,y,dy);
      }
      terror = fabs(tc0 - y[TC]) + fabs(te0 - y[TE]);
      printf("              final Tc = %.1f(K), Te = %.1f(K)\n", y[TC], y[TE]);
      printf("   temperature error (del(Tc) + del(Te)) = %.3f(K)\n", terror);
      fprintf(printfile, "              final Tc = %.1f(K), Te = %.1f(K)\n",
               y[TC], y[TE]);
      fprintf(printfile, "   temperature error (del(Tc) + del(Te)) = %.3f(K)\n",
               terror);
      iter++;
   } while((terror >= epsilon) && (iter <= maxit));
   if(iter > maxit) {
      printf("\n no convergence within %d iterations\n", maxit);
      fprintf(printfile,"\n no convergence within %d iterations\n", maxit);
   }
/*     *** setup output matrices **************************************/
   theta = 0.0;
   y[QK] = 0.0;
   y[QR] = 0.0;
   y[QH] = 0.0;
   y[WC] = 0.0;
   y[WE] = 0.0;
   y[W] = 0.0;

   filmat(0,y,dy,var,dvar);
   for(i=1;i<COL;i++) {
      for(j=0;j<step;j++)
         rk4(dadiab,7,&theta,dtheta,y,dy);
      filmat(i,y,dy,var,dvar);
   }
}