
// cpp.c  -  Modification from the original file: cpropep.c            //
// from SourceForge Project:  Rocket Workbench Project                 //
// http://sourceforge.net/projects/rocketworkbench/                    //

// cpropep.c  -  Calculation of Complex Chemical Equilibrium           //
// $Id: cpropep.c,v 1.2 2000/07/19 02:13:03 antoine Exp $              //
// Copyright (C) 2000                                                  //
//    Antoine Lefebvre <antoine.lefebvre@polymtl.ca>                   //
//    Mark Pinese <pinese@cyberwizards.com.au>                         //
//                                                                     //
// Licensed under the GPLv2                                            //

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <malloc.h>
#include <time.h>

#include "getopt.h"
#include "load.h"
#include "equilibrium.h"
#include "performance.h"
#include "derivative.h"
#include "thermo.h"

#include "print.h"

#include "conversion.h"
#include "compat.h"
#include "return.h"

#define version "1.0"
#define date    "10/07/2000"

#define CHAMBER_MSG     "Time spent for computing chamber equilibrium"
#define FROZEN_MSG      "Time spent for computing frozen performance"
#define EQUILIBRIUM_MSG "Time spent for computing equilibrium performance"

#define CONF_FILE "cpropep.conf"

#define TIME(function, msg) timer = clock(); function;\
                            fprintf(outputfile,\
                                    "%s: %f s\n\n", msg,\
                                    (double)(clock() - timer)/CLOCKS_PER_SEC)


//#undef TIME
//#define TIME(function, msg) function;

#define MAX_CASE 10

typedef enum _p
{
  SIMPLE_EQUILIBRIUM,
  FIND_FLAME_TEMPERATURE,
  FROZEN_PERFORMANCE,
  EQUILIBRIUM_PERFORMANCE
} p_type;

char case_name[][80] = {
  "Fixed pressure-temperature equilibrium",
  "Fixed enthalpy-pressure equilibrium - adiabatic flame temperature",
  "Frozen equilibrium performance evaluation",
  "Shifting equilibrium performance evaluation"
};

char thermo_file[FILENAME_MAX] = "thermo.dat";
char propellant_file[FILENAME_MAX] = "propellant.dat";


typedef struct _case_t
{
  p_type p;

  bool temperature_set;
  bool pressure_set;
  bool exit_condition_set;

  double           temperature;
  double           pressure;
  exit_condition_t exit_cond_type;
  double           exit_condition;
  
} case_t;


// CPROPEP DLL Exported Functions

#define EXPORT __declspec(dllexport) __stdcall

//BOOL APIENTRY DllMain(HANDLE hModule, DWORD  ul_reason_for_call, LPVOID lpReserved){return TRUE;}

typedef struct _CComp{
	double Mass;
    int    Code;
    int    H;
	} CComp;
  
typedef struct _CPro{
    // output data
	equilib_prop_t     Prop[3];
	performance_prop_t Perf[3];
	double             Density;
	// input data
    double ChamberTemperature;
    double ChamberPressure;
    double ExitPressure;
    int   Mode;
	int   ExitMode;
    int   Components;
    int   H;
	CComp Comp[MAX_COMP];
    } CPro;

void EXPORT InitCP(char *propellant_file, char *thermo_file){
	
	load_thermo(thermo_file);
	load_propellant(propellant_file);
	
	}

void EXPORT CleanUpCP(){
	
	free(propellant_list);
	free(thermo_list);
	
	}

void EXPORT NameCP(int c, char *t){
	
	if (propellant_list == NULL || t == NULL || c < 0 || c >= num_propellant)
		return;
	strcpy(t, (propellant_list + c)->name);
	
	}

void EXPORT CompCP(int c, char *t){
	
	int n;
	if (propellant_list == NULL || t == NULL || c < 0 || c >= num_propellant)
		return;
	for (n = 0; n < 6; n++){
		if ((propellant_list + c)->coef[n] != 0)
			t+=sprintf(t, "%d%s ", (propellant_list + c)->coef[n], symb[(propellant_list + c)->elem[n]]);
		}
	
	}

int EXPORT CodeCP(char *t){
	
	int c, n;
	n = strlen(t);
	if (propellant_list == NULL || t == NULL)
		return -1;
	for (c = 0; c < num_propellant; c++)
		if (!strncmp(t, (propellant_list + c)->name, n))
			return c;
	return -1;
	
	}

void EXPORT ExecuteCP(CPro *c){
	
	int i;
	equilibrium_t *equil, *frozen, *shifting;
	
	if (propellant_list == NULL || thermo_list == NULL || c == NULL)
		return;
	
	outputfile = stdout; //NULL
	errorfile = stderr; //NULL
	
	equil = (equilibrium_t *) malloc(sizeof (equilibrium_t));
	initialize_equilibrium(equil);
	frozen   = (equilibrium_t *) malloc(sizeof(equilibrium_t) * 3);
	shifting = (equilibrium_t *) malloc(sizeof(equilibrium_t) * 3);
	for (i = 0; i < 3; i++){
		initialize_equilibrium(frozen + i);
		initialize_equilibrium(shifting + i);
		}
	
	for (i = 0; i < c->Components; i++)
		add_in_propellant(equil, c->Comp[i].Code, GRAM_TO_MOL(c->Comp[i].Mass, c->Comp[i].Code));
	
	compute_density(&(equil->propellant));
	list_element(equil);
	list_product(equil);
	
	switch (c->Mode){
		
		case SIMPLE_EQUILIBRIUM:
		
			equil->properties.T = c->ChamberTemperature;
			equil->properties.P = c->ChamberPressure;
			//print_propellant_composition(equil);
			c->Density = equil->propellant.density;
			//resolve equilibrium
			if (equilibrium(equil, TP)) break;
			//print_product_properties(equil, 1);
			//print_product_composition(equil, 1);
			break;
		
		case FIND_FLAME_TEMPERATURE:
		
			equil->properties.P = c->ChamberPressure;
			//print_propellant_composition(equil);
			c->Density = equil->propellant.density;
			//resolve equilibrium
			if (equilibrium(equil, HP)) break;
			//print_product_properties(equil, 1);
			//print_product_composition(equil, 1);
			break;
		
		case FROZEN_PERFORMANCE:
		
			equil->properties.T = c->ChamberTemperature;
			equil->properties.P = c->ChamberPressure;
			copy_equilibrium(frozen, equil);
			//print_propellant_composition(frozen);
			c->Density = frozen->propellant.density;
			//resolve equilibrium
			if (equilibrium(frozen, HP)) break;
			//execute frozen performance
			frozen_performance(frozen, (exit_condition_t) c->ExitMode, c->ExitPressure);
			//print_product_properties(frozen, 3);
			c->Prop[0] = (frozen)->properties;
			c->Prop[1] = (frozen + 1)->properties;
			c->Prop[2] = (frozen + 2)->properties;
			//print_performance_information(frozen, 3);
			c->Perf[0] = (frozen)->performance;
			c->Perf[1] = (frozen + 1)->performance;
			c->Perf[2] = (frozen + 2)->performance;
			//print_product_composition(frozen, 3);
			break;
		
		case EQUILIBRIUM_PERFORMANCE:
		
			equil->properties.T = c->ChamberTemperature;
			equil->properties.P = c->ChamberPressure;
			copy_equilibrium(shifting, equil);
			//print_propellant_composition(shifting);
			c->Density = shifting->propellant.density;
			//resolve equilibrium
			if (equilibrium(shifting, HP)) break;
			//execute shifting performance
			shifting_performance(shifting, (exit_condition_t) c->ExitMode, c->ExitPressure);
			//print_product_properties(shifting, 3);
			c->Prop[0] = (shifting)->properties;
			c->Prop[1] = (shifting + 1)->properties;
			c->Prop[2] = (shifting + 2)->properties;
			//print_performance_information(shifting, 3);
			c->Perf[0] = (shifting)->performance;
			c->Perf[1] = (shifting + 1)->performance;
			c->Perf[2] = (shifting + 2)->performance;
			//print_product_composition(shifting, 3);
			break;
		
		}
	
	free (equil);
	free (frozen);
	free (shifting);
	
	}
