/***************************************************************************/
/*                                                                         */
/* The sor06.C main program                                                */
/*  Version: 1.0                                                           */
/*  Date: 11 Aug. 2018                                                     */
/*                                                                         */
/*  Primary Authors: M.S. Olufsen                                          */
/*  Key Contributers: M.U. Qureshi & M.J. Colebank                         */
/*  Institute: North Carolina State University, Raleigh, NC, USA           */
/*  Funded by: NSF-DMS # 1615820                                           */
/*                                                                         */
/***************************************************************************/

// $Id: sor06.C,v 1.0 08/11/2018 16:00:42 MUQ $

#include "sor06.h"
#include "tools.h"
#include "arteries.h"
#include <stdlib.h> // for atoi, atof

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cstring>

using namespace std;

// The vessel-network is specified and initialized, and the flow and
// pressures are to be determed. All global constants must be defined
// in the header file. That is the number of dots per cm, and the number
// of vessels.
int main(int argc, char *argv[])
{
  double tstart, tend, finaltime;

    
    double f1   = 0;
    double f2   = 1;
    
    double R41, R61, R81, R101, R111, R121, R141, R161, R181, R191, R201;
    double R42, R62, R82, R102, R112, R122, R142, R162, R182, R192, R202;
    double C4, C6, C8, C10, C11, C12, C14, C16, C18, C19, C20;
    
    R41 = 510.4072231; R61 = 388.0002164; R81 = 538.8259481; R101 = 466.9618518; R111= 1807.757766; R121 = 1173.676967;
    R141 = 266.0796751; R161 = 1041.232548; R181 = 107.1125395; R191 = 116.974148; R201 = 83.43084455;
    
    R42 = 2041.628892; R62 = 1552.000866; R82 = 2155.303792; R102 = 1867.847407; R112 = 7231.031065; R122 = 4694.70787;
    R142 = 1064.3187;R162= 4164.930191; R182 = 428.4501579; R192 = 467.896592; R202 = 333.7233782;
    
    C4 = 0.005254628; C6 = 0.006912367; C8 = 0.004977489; C10 = 0.00574351; C11 = 0.001483606; C12 = 0.002389967;
    C14 = 0.01054214; C16 = 0.00269397; C18 = 0.026187868; C19 = 0.023980077; C20 = 0.033621248;
    
    double f3, r1, r2, c1;
    int HB, cycles, id;
    
    
    if (argc != 8) //argv[0] is the name of the program, here sor06
    {
        printf("Not enough input arguments, noargc %d and they are %s\n", argc, argv[0]);
        return 1;
    }
    
    f3 = atof(argv[1]);
    r1 = atof(argv[2]);
    r2 = atof(argv[3]);
    c1 = atof(argv[4]);
    HB = atoi(argv[5]);
    cycles = atoi(argv[6]);
    id = atoi(argv[7]);

	char namepu1 [20];
    char namepu2 [20];
    char namepu3 [20];
	
    sprintf(namepu1, "pu1_%d.2d", id);
    sprintf(namepu2, "pu2_%d.2d", id);
    sprintf(namepu3, "pu3_%d.2d", id);
    
	
    FILE *fpu1 = fopen (namepu1, "w");
    FILE *fpu2 = fopen (namepu2, "w");
    FILE *fpu3 = fopen (namepu3, "w");
		
  // Workspace used by bound_bif
  for(int i=0; i<18; i++) fjac[i] = new double[18];

  clock_t cc1 = clock();        // Only used when timing the program.
  nbrves    = 21;             // The number of vessels in the network.
  tstart    = 0.0;            // Starting time.
  finaltime = HB*Period;       // Final end-time during a simulation. numHeartBeats = 6
  tend      = (HB-cycles)*Period;       // Timestep before the first plot-point is reached.
    

  Tube   *Arteries[nbrves];                    // Array of blood vessels.
    
    
    // Initialization of the Arteries.

   
    // Definition of Class Tube: (Length, topradius, botradius, *LeftDaughter, *RightDaughter, points, init, K, f1, f2, f3, R1, R2,  CT);

    
    // =================ORIGINAL NETWORK Control (Poiseuille based flow distribution for WK)==============================
   
    Arteries[20] = new Tube( 0.055, 0.018, 0.018, 0, 0, 40,0,0,f1,f2,f3, r1*R201, r2*R202, c1*C20);//RH-PW //G5
    Arteries[19] = new Tube( 0.178, 0.022, 0.022, 0, 0, 40,0,0,f1,f2,f3, r1*R191, r2*R192, c1*C19);//RH-PW //G5
    Arteries[18] = new Tube( 0.177, 0.015, 0.015, 0, 0, 40,0,0,f1,f2,f3, r1*R181, r2*R182, c1*C18);//RH-PW //G4
    
    Arteries[17] = new Tube( 0.469, 0.024, 0.024, Arteries[ 19], Arteries[ 20], 40,0,0,f1,f2,f3, 0, 0, 0);
    
    Arteries[16] = new Tube( 0.302, 0.015, 0.015, 0, 0, 40,0,0,f1,f2,f3, r1*R161, r2*R162, c1*C16);//RH-PW //G3
    
    Arteries[15] = new Tube( 0.083, 0.025, 0.025, Arteries[ 17], Arteries[ 18], 40,0,0,f1,f2,f3, 0, 0, 0);
    
    Arteries[14] = new Tube( 0.184, 0.019, 0.019, 0, 0, 40,0,0,f1,f2,f3, r1*R141, r2*R142, c1*C14);//RH-PW //G2
    
    Arteries[13] = new Tube( 0.081, 0.026, 0.026, Arteries[ 15], Arteries[ 16], 40,0,0,f1,f2,f3, 0, 0, 0);
    
    Arteries[12] = new Tube( 0.062, 0.014, 0.014, 0, 0, 40,0,0,f1,f2,f3, r1*R121, r2*R122, c1*C12);//LH-PW //G4
    
    Arteries[11] = new Tube( 0.140, 0.015, 0.015, 0, 0, 40,0,0,f1,f2,f3, r1*R111, r2*R112, c1*C11);//LH-PW //G4
    
    Arteries[10] = new Tube( 0.069, 0.016, 0.016, 0, 0, 40,0,0,f1,f2,f3, r1*R101, r2*R102, c1*C10);//LH-PW //G3
    
    Arteries[9] = new Tube( 0.262, 0.020, 0.020, Arteries[ 11], Arteries[ 12], 40,0,0,f1,f2,f3, 0, 0, 0);
    
    Arteries[8] = new Tube( 0.177, 0.017, 0.017, 0, 0, 40,0,0,f1,f2,f3, r1*R81, r2*R82, c1*C8);//LH-PW //G2
    
    Arteries[7] = new Tube( 0.311, 0.023, 0.023, Arteries[ 9], Arteries[ 10], 40,0,0,f1,f2,f3, 0, 0, 0);
    
    Arteries[6] = new Tube( 0.212, 0.017, 0.017, 0, 0, 40,0,0,f1,f2,f3, r1*R61, r2*R62, c1*C6);//RH-PW //G1
    
    Arteries[5] = new Tube( 0.202, 0.032, 0.032, Arteries[ 13], Arteries[ 14], 40,0,0,f1,f2,f3, 0, 0, 0);
    
    Arteries[4] = new Tube( 0.052, 0.013, 0.013, 0, 0, 40,0,0,f1,f2,f3, r1*R41, r2*R42, c1*C4);//LH-PW //G1
    
    Arteries[3] = new Tube( 0.241, 0.024, 0.024, Arteries[ 7], Arteries[ 8], 40, 0, 0,f1,f2,f3, 0, 0, 0);
    
    Arteries[2] = new Tube( 0.372, 0.037, 0.037, Arteries[ 5], Arteries[ 6], 40, 0, 0,f1,f2,f3, 0, 0, 0);
    
    Arteries[1] = new Tube( 0.445, 0.026, 0.026, Arteries[ 3], Arteries[ 4], 40, 0, 0,f1,f2,f3, 0, 0, 0);
    
    Arteries[0] = new Tube( 0.410, 0.047, 0.047, Arteries[ 1], Arteries[ 2], 40, 1, 0,f1,f2,f3, 0, 0, 0);

    
    

  // In the next three statements the simulations are performed until
  // tstart = tend. That is this is without making any output during this
  // first period of time. If one needs output during this period, these three
  // lines should be commented out, and the entire simulation done within the
  // forthcomming while-loop.

  // Solves the equations until time equals tend.
  solver (Arteries, tstart, tend, k);
  tstart = tend;
  tend = tend + Deltat;

    
  // The loop is continued until the final time
  // is reached. If one wants to make a plot of
  // the solution versus x, tend is set to final-
  // time in the above declaration.
  while (tend <= finaltime)
  {
    for (int j=0; j<nbrves; j++)
    {
      int ArtjN = Arteries[j]->N;
      for (int i=0; i<ArtjN; i++)
      {
        Arteries[j]->Qprv[i+1] = Arteries[j]->Qnew[i+1];
        Arteries[j]->Aprv[i+1] = Arteries[j]->Anew[i+1];
      }
    }

    // Solves the equations until time equals tend.
    solver (Arteries, tstart, tend, k);
	
	
    // A 2D plot of P(x_fixed,t) is made. The resulting 2D graph is due to
    // the while loop, since for each call of printPt only one point is set.
      Arteries[ 0] -> printALLxt (fpu1, tend, 0);
      Arteries[ 1] -> printALLxt (fpu2, tend, 0);
      Arteries[ 2] -> printALLxt (fpu3, tend, 0);
      
     
      
    // The time within each print is increased.
    tstart = tend;
    tend   = tend + Deltat; // The current ending time is increased by Deltat.
  }
  
    
  // The following statements is used when timing the simulation.
  fprintf(stdout,"nbrves = %d, Lax, ", nbrves);
  clock_t cc2 = clock(); // FIXME clock() may wrap after about 72 min.
  int tsec = (int) ((double) (cc2-cc1)/CLOCKS_PER_SEC + 0.5);
  fprintf(stdout,"cpu-time %d:%02d\n", tsec / 60, tsec % 60);
  
  

  // In order to terminate the program correctly the vessel network and hence
  // all the vessels and their workspace are deleted.
  for (int i=0; i<nbrves; i++) delete Arteries[i];

  // Matrices and arrays are deleted
  for (int i=0; i<18; i++) delete[] fjac[i];
    fclose (fpu1);
    fclose (fpu2);
    fclose (fpu3);

  return 0;
}
