/******************************************************************
 * FUNCTION: simple -Stirling machine analysis with non-ideal     *
 *         heat exchanger and appendix losses                     *
 * FILENAME simple.c                                              *
 * PROGRAMMER I.Urieli    (FORTRAN)                               *
 *            E. Malroy   (C Translation)                         *
 * DATE:      12/01/94  (FORTRAN), 06/16/96 (C Translation)       *
 * LAST MODIFICATION: 09/05/96 by I Urieli                        *
 * INCLUDE:  <stdio.h>, <math.h>, "adiabatic.h", "define.h",      *
 *    "simple.h"                                                  *
 * GLOBAL VARIABLES:                                              *
 *    freq   - frequency of engine                                *
 *    tk     - cooler temperature (K)                             *
 *    th     - heater temperature (K)                             *
 *    tr     - regenerator temperature                            *
 *    cqwr   - regenerator wall thermal conductance (Watts/K)     *
 * PROTOTYPES:                                                    *
 *    void simple(void)                                           *
 * PRE:  freq, tk, th, tr, cqwr                                   *
 * POST:  th, tk, qwrl, qrloss, shutls, enthls, sealls, dwork     *
 ******************************************************************/
#include <stdio.h>
#include <math.h>
#include "adiabatic.h"
#include "define.h"
#include "simple.h"

void simple(void)

{

double  epsilon = 1.0; /* Temperature error (deg. K)        */
double var[ROWV][COL],dvar[ROWD][COL];
/* var,dvar: matrices of variables and derivatives          */

double twh;      /* heater wall temperature (K)             */
double tgh;      /* heater average gas temperature (K)      */
double twk;      /* kooler wall temperature (K)             */
double tgk;      /* kooler average gas temperature (K)      */
double qrloss;   /* regen net enthalpy loss (J/rad)         */
double qwrl;     /* regen wall heat leakage (J/rad)         */
double dwork;  /* pressure drop available work loss (J/rad) */
double shutls; /* shuttle loss (Watts)                      */
double enthls; /* gap enthalpy (pumping) loss (Watts)       */
double sealls; /* displacer seal PV loss (Watts)            */

   printf("\n\n The \'Simple\' Stirling cycle analysis...");
   fprintf(printfile, "\n The \'Simple\' Stirling cycle analysis...");
   twk = tk;
   twh = th;
   do
   {
      adiab(var,dvar);
      hotsim(var,twh,&tgh);
      kolsim(var,twk,&tgk);
      if((fabs(th - tgh) + fabs(tk - tgk)) < epsilon)  break;
      th = tgh;
      tk = tgk;
      tr = (th - tk)/log(th/tk);
   }
   while(1);
   printf("\nHeat Exchanger analysis converged...\n");
   printf(" Ideal Adiabatic output: Th = %.4f, Tk = %.4f (K)\n",th,tk);

   fprintf(printfile,"\nHeat Exchanger analysis converged...\n");
   fprintf(printfile,"Ideal Adiabatic output: Th = %.4f, Tk = %.4f (K)\n",
             th,tk);
   prntad(var,dvar);

   printf("\n Regenerator simple analysis...\n");
   printf(" ideal heat accepted(Watts) %.3f\n",var[QH][36]*freq);
   printf(" ideal indicated power(Watts) %.3f\n",var[W][36]*freq);

   fprintf(printfile,"\n Regenerator simple analysis...\n");
   fprintf(printfile," ideal heat accepted(Watts) %.3f\n",var[QH][36]*freq);
   fprintf(printfile," ideal indicated power(Watts) %.3f\n",var[W][36]*freq);

   regsim(var,&qrloss);
   if(th < tk) {    /* if refrigerator... */
      qrloss = -qrloss;
   }
   printf("regenerator net enthalpy loss(Watts) %.3f\n", qrloss*freq);
   fprintf(printfile,"regenerator net enthalpy loss(Watts) %.3f\n",
             qrloss*freq);
   qwrl = cqwr*(twh - twk)/freq;
   printf("regenerator wall heat leakage(Watts) %.3f\n", qwrl*freq);
   fprintf(printfile,"regenerator wall heat leakage(Watts) %.3f\n", qwrl*freq);

   wrksim(var,dvar,&dwork);

   printf("pressure drop available work loss(Watts) %.3f\n",dwork*freq);
   fprintf(printfile,"pressure drop available work loss(Watts) %.3f\n",
       dwork*freq);

 /***Appendix gap losses - relevant if displacer or hot piston data available */
/*   aploss(var,shutls,enthls,sealls)
   printf("  shuttle loss (Watts)  %.3f\n", shutls);
   printf("  gap entalpy (pumping) loss (Watts)  %.3f\n", enthls);
   printf("  displacer seal PV loss (Watts)  %.3f\n", sealls);
   fprintf(printfile,"  shuttle loss (Watts)  %.3f\n", shutls);
   fprintf(printfile,"  gap entalpy (pumping) loss (Watts)  %.3f\n", enthls);
   fprintf(printfile,"  displacer seal PV loss (Watts)  %.3f\n", sealls);  */

}