/******************************************************************
* FUNCTION: dadiab - evaluate Ideal Adabaatic Model derivatives   *
* FILENAME dadiab.c                                               *
* PROGRAMMER I.Urieli    (FORTRAN)                                *
*            E. Malroy   (C Translation)                          *
* DATE:      11/30/94  (FORTRAN), 03/12/96 (C Translation)        *
* LAST MODIFIED: 08/17/96 by I.Urieli                             *
* INCLUDE:  <stdio.h>,<math.h>, "adiabatic.h","define.h"  *
* GLOBAL VARIABLES: none                                          *
* PROTOTYPES:                                                     *
* PRE:   none                                                     *
* POST:  none                                                     *
******************************************************************/
#include <stdio.h>
#include <math.h>
#include "define.h"
#include "adiabatic.h"

void dadiab(double theta,   /* crank angle (radians)             */
            double y[ROWV], /* vector of dependant variables     */
            double dy[ROWD]) /* vector of derivatives            */
{
 /* variables */
   double vot, dpop, hold;
   int i;
/*     ***volumes and volume derivatives*****************************/
   volume(theta,&y[VC],&y[VE],&dy[VC],&dy[VE]);
/*     ***pressure and pressure derivative***************************/
   vot = vk/tk + vr/tr + vh/th;
   y[P] = mgas*rgas/(y[VC]/y[TC] + vot + y[VE]/y[TE]);
   hold = (y[VC]/(y[TCK]*gama)+vot+y[VE]/(y[THE]*gama));
   dy[P] = -y[P]*(dy[VC]/y[TCK]+dy[VE]/y[THE])/hold;
   dpop = dy[P]/y[P];
/*     ***masses and mass accumulations******************************/
   y[MC] = y[P]*y[VC]/(rgas*y[TC]);
   y[MK] = y[P]*vk/(rgas*tk);
   y[MR] = y[P]*vr/(rgas*tr);
   y[MH] = y[P]*vh/(rgas*th);
   y[ME] = y[P]*y[VE]/(rgas*y[TE]);
   dy[MC] = (y[P]*dy[VC] + y[VC]*dy[P]/gama)/(rgas*y[TCK]);
   dy[ME] = (y[P]*dy[VE] + y[VE]*dy[P]/gama)/(rgas*y[THE]);
   dy[MK] = y[MK]*dpop;
   dy[MR] = y[MR]*dpop;
   dy[MH] = y[MH]*dpop;
/*     ***mass flows*************************************************/
   y[GACK] = -dy[MC];
   y[GAKR] = y[GACK] - dy[MK];
   y[GAHE] = dy[ME];
   y[GARH] = y[GAHE] + dy[MH];
/*     ***conditional temperatures************************************/
   y[TCK] = tk;
   if(y[GACK] >= 0.0)
      y[TCK] = y[TC];
   y[THE] = y[TE];
   if(y[GAHE] >= 0.0)
      y[THE] = th;
/*     ***working space temperature derivatives**********************/
   dy[TC] = y[TC]*(dpop + dy[VC]/y[VC] - dy[MC]/y[MC]);
   dy[TE] = y[TE]*(dpop + dy[VE]/y[VE] - dy[ME]/y[ME]);
/*     ***energy*****************************************************/
   dy[QK] = vk*dy[P]*cv/rgas - cp*(y[TCK]*y[GACK] - tk*y[GAKR]);
   dy[QR] = vr*dy[P]*cv/rgas - cp*(tk*y[GAKR] - th*y[GARH]);
   dy[QH] = vh*dy[P]*cv/rgas - cp*(th*y[GARH] - y[THE]*y[GAHE]);
   dy[WC] = y[P]*dy[VC];
   dy[WE] = y[P]*dy[VE];
   dy[W] = dy[WC] + dy[WE];
   y[W] = y[WC] + y[WE];

}