/************************************************************
 * FUNCTION: regsim -Evaluate the effectiveness and         *
 *   performance of the regenerator                         *
 *                                                          *
 * FILENAME regsim.c                                        *
 * PROGRAMMER I.Urieli    (FORTRAN)                         *
 *            E. Malroy   (C Translation)                   *
 * DATE:      12/01/94  (FORTRAN), 06/15/96 (C Translation) *
 * LAST MODIFICATION:  09/05/96 by I.Urieli                 *
 * INCLUDE:  <stdio.h>,<math.h>,"adiabatic.h","define.h"    *
 *           "simple.h"                                     *
 * GLOBAL VARIABLES:                                        *
 *    omega: operating frequency (rads/sec)                 *
 *    tr: effective regenerator temperature (K)             *
 *    dr: regenerator hydraulic diameter (m)                *
 *    awgr: regenerator internal wetted area (sq.m)         *
 *    ar: regenerator internal free flow area (sq.m         *
 *    rmatrx:  char that tells regenerator type ("f" - foil *
 *             matrix and "m" - mesh matrix)                *
 *    printfile:  points to output file                     *
 * PROTOTYPES:                                              *
 *    void regsim(double [][], double *)                    *
 * PRE:   omega, tr, dr, awgr, ar, printfile                *
 * POST:  var[gahe][37], qrloss                             *
 ************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "adiabatic.h"
#include "define.h"
#include "simple.h"

void regsim(double var[ROWV][COL],  /* matrix of output variables */
            double *pt_qrloss)      /* Amount of heat that must be
               added to the heater due to the non-ideal regenerator*/
{
  double gar[37];  /* mass flow through regenerator (kg/s) */
  double gr;       /* mass flux through regenerator (kg/(m**2*s)  */
  double re[37];   /* Reynolds number                             */
  double sumre;    /* sumation of Reynolds number                 */
  double reavg;    /* average Reynolds number over the cycle      */
  double remax;    /* maximum Reynolds number over the cycle      */
  double st;       /* Stanton number                              */
  double ht;       /* heat transfer coefficient (W/(m^2*K))       */
  double effect;   /* regenerator effectiveness                   */
  double qreg[37]; /* regenerator heat flow vector (used with MinMax) */
  double fr;       /* Reynolds friction factor                    */
  double mu;       /* dynamic viscosity (kg/(m*s))                */
  double kgas;     /* gas thermal conductivity (W/(m**2*K))       */
  double ntu;      /* Number of transfer units                    */
  double qrmin;    /* minimum amount of heat added due to         */
                   /* non-ideal regenerator over the cycle        */
  double qrmax;    /* maximum amount of heat added due to         */
                   /* non-ideal regenerator over the cycle        */
  double qrloss;   /* Amount of heat added to heater due to       */
                   /* non-ideal regenerator                       */
  int i;           /* index for "for" loop                        */

/**** Reynolds number over the cycle ****/
   for(i=0; i<37; i++) {
      gar[i] = (var[GAKR][i] + var[GARH][i])*omega/2.0;
      gr = gar[i]/ar;
      reynum(tr,gr,dr,&mu,&kgas,&re[i]);
   }

/**** average and maximum Reynolds number****/
   sumre = 0.0;
   remax = re[0];
   for( i=0; i<36; i++) {
      sumre = sumre + re[i];
      if(re[i] > remax)
         remax = re[i];
   }
   reavg = sumre/36.;

/**** Stanton number, number of transfer units and regen effectiveness ****/
   if(rmatrx == 'm')
      matfr(reavg,&st,&fr);
   else if(rmatrx == 'f')
      foilfr(dr,mu,reavg,&st,&ht,&fr);
   else {
      printf("\nundefined regenerator matrix\n");
      fprintf(printfile,"\nundefined regenerator matrix\n");
      exit(0);
   }
   ntu = st*awgr/(2.0*ar);
   effect = ntu/(ntu+1);

/**** calculate Qrloss ****/
   for(i=0; i<37; i++)
      qreg[i] = var[QR][i];
   minmax(qreg,36,&qrmin,&qrmax);

/**** check in-out qrmin, qrmax ****/
   qrloss = (1.0 - effect)*(qrmax-qrmin);
   *pt_qrloss = qrloss;

/**** output ****/
   printf(" average Reynolds number %.3f\n",reavg);
   printf(" maximum Reynolds number %10.3lf\n",remax);
   printf(" Stanton number(Average Re) %.3f\n",st);
   printf(" Number of transfer units %.3f\n",ntu);
   printf(" regenerator effectivenessct %.3f\n",effect);
   printf(" Qrloss(Joules) %.4f\n",qrloss);

   fprintf(printfile," average Reynolds number %.3f\n",reavg);
   fprintf(printfile," maximum Reynolds number %.3f\n",remax);
   fprintf(printfile," Stanton number(Average Re) %.3f\n",st);
   fprintf(printfile," Number of transfer units %.3f\n",ntu);
   fprintf(printfile," regenerator effectivenessct %.3f\n",effect);
   fprintf(printfile," Qrloss(Joules) %.4f\n",qrloss);
}

/******************************************************************
* FUNCTION: minmax -find both ymin and ymax values in a real array*
*                   Y of n values                                 *
*                                                                 *
* FILENAME regsim.c                                               *
* PROGRAMMER I.Urieli    (FORTRAN)                                *
*            E. Malroy   (C Translation)                          *
* DATE:      12/01/94  (FORTRAN), 06/15/96 (C Translation)        *
* LAST MODIFICATION:  08/21/96 Eric Malroy                        *
* INCLUDE:  <stdio.h>,<math.h>,"regsim.h"                         *
* GLOBAL VARIABLES: None                                          *
* PROTOTYPES:                                                     *
*    void minmax(double [], int, double *, double *)              *
* PRE:   None                                                     *
* POST:  *ymin, *ymax                                             *
*******************************************************************/
void minmax(double y[],  /* array to find min and max             */
            int n,      /* int that tells the number in the array */
            double *pt_ymin, /* the minimum in the array          */
            double *pt_ymax) /* the maximum in the array          */
{

 double y_min, y_max;  /* minimum or maximum y value */
 int i;

   y_min = y[0];
   y_max = y[0];
   for( i = 1; i<= n; i++)
   {
      if(y[i] < y_min)
         y_min = y[i];
      else if (y[i] > y_max)
         y_max = y[i];
   }
   *pt_ymin = y_min;
   *pt_ymax = y_max;
}