#ifndef lint
static char *rcs="$Header: /mit/hst/src/sim/RCS/mktbl.c,v 2.1 90/08/18 18:40:26 tldavis Exp $";
#endif
/*
 *	$Source: /mit/hst/src/sim/RCS/mktbl.c,v $
 *
 *  This program gemerates capacitance file and dc/dt file based on model
 *  first suggested by Bob Sah in his original notes on the simulator.  This
 *  fixes a bug which occurred when the original elastance minimum was 0.
 *  This version has minimum elastance of 0.1, maximum of 1, as in paper
 *  in CRC Biomedical Engineering  ( I'll find the reference later.)
 *  Relaxation occurs in about 1/2 the time as contraction, both similar
 *  to sinusoidal elastance.
 *  Careful to round to nearest integer!  Using floor(x+ 0.5) to do this.
 *  dcap is hand-done derivative of cap with respect to i.
 */

#include <stdio.h>
#include <math.h>
#include "CVDefs.h"

main()
{
  register double elast;
  register double cap, dcap;
  register int i;
  double pi = 3.14159265;
  FILE *tblc, *tbldcdt;
  
  tblc = fopen("tblc", "w+");  tbldcdt = fopen("tbldcdt", "w+");
  if (!tblc || !tbldcdt) {
    printf("Couldn't open file to write!\n");
    exit(1);
  }
  fprintf(tblc, "%d,\n\n\n", 90000);    /* these are the max value for c */
  fprintf(tbldcdt, "%d,\n\n\n", 90000);

  for (i=0; i < SYSTOLESEGS; i++) { /* during slow contraction */
    elast = .1 + .9*(cos((i - SYSTOLESEGS)*pi/SYSTOLESEGS) + 1)/2;
    cap = 10000.0*(1.0/elast - 1.0);
    fprintf(tblc, "%d,\n", (int) floor(cap + 0.5));

    dcap = 10000.0*(pi/SYSTOLESEGS)*0.45 * 
      sin((i-SYSTOLESEGS)*pi/SYSTOLESEGS) / (elast*elast);
    fprintf(tbldcdt, "%d, \n", (int) floor(dcap + 0.5));
  }

  for (i=SYSTOLESEGS; i < TBLSEGSPLUS2; i++) { /* during fast relaxation */
    elast = .1 + .9*(cos((i - SYSTOLESEGS)*pi*2/SYSTOLESEGS) + 1)/2;
    cap = 10000.0*(1.0/elast - 1.0);
    fprintf(tblc, "%d,\n", (int) floor(cap + 0.5));

    dcap = 10000.0*(2*pi/SYSTOLESEGS)*0.45 * 
      sin((i-SYSTOLESEGS)*2*pi/SYSTOLESEGS) / (elast*elast);
    fprintf(tbldcdt, "%d,\n", (int) floor(dcap + 0.5));
  }
  exit (0);
}