// The Aries Systems Code, SiC module
// LCC cleaned up formatting 8/2010

#include <cstdlib>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string.h>
#include <math.h>

using std::string;
using namespace std;

#include "Part.h"
#include "DesignPoint.h"
#include "CostingAccount.h"
#include "Geometry.h"
#include "Cost.h"

// Decplare global variable available everywhere
string blanketmode, dispmode, manifold;

void system_generator(double (*chuck)[100], int ndp, 				string dispmode, 
					  string blanketmode, 	string plasmaswitch, 	string manifold);

extern"C"
{
void sys_(
	double *rq,		double *aspq,	double *btq,	double *aipq,		double *betanq,
    double *qmhdq,	double *qcylq,	double *deltaq,	double *anq,		double *atq,
    double *fracgwq,double *Qq,		double *akq,	double *hfactorq,	double *taujq,
    double *taueq,	double *tpoteq,	double *xneq,	double *tempaveq,	double *tflatq,   
    double *pltohq,	double *xnwallq,double *pbremq,	double *fbsq,		double *etacdq,
    double *powcdq,	double *pauxq,	double *pcyclq,	double *pohmq,		double *plineq,
    double *betafq,	double *zeffq,	double *tratq,	double *dratq,		double *powcd1q, 
    double *powcd2q,double *powcd3q,double *powcd4q,double *fimp1q,		double *fimp2q,  
    double *fimp3q,	double *temp0q,	double *fcdq,	double *fniq,		double *fheq,
    double *fdtq,	double *wtotq,	double *dfcsq,	double *vloopq,		double *dptavq,
    double *tptavq,	double *betapq,	double *betatq,	double *palphaq,	double *ploplhq);
     
void sysengr_(
	double *rq,		double *aspq,	double *btq,	double *aipq,		double *betanq,
    double *qmhdq,	double *qcylq,	double *deltaq,	double *anq,		double *atq,
    double *fracgwq,double *Qq,		double *akq,	double *hfactorq,	double *taujq,
    double *taueq,	double *tpoteq,	double *xneq,	double *tempaveq,	double *tflatq,   
    double *pltohq,	double *xnwallq,double *pbremq,	double *fbsq,		double *etacdq,
    double *powcdq,	double *pauxq,	double *pcyclq,	double *pohmq,		double *plineq,
    double *betafq,	double *zeffq,	double *tratq,	double *dratq,		double *powcd1q, 
    double *powcd2q,double *powcd3q,double *powcd4q,double *fimp1q,		double *fimp2q,  
    double *fimp3q,	double *temp0q,	double *fcdq,	double *fniq,		double *fheq,
    double *fdtq,	double *wtotq,	double *dfcsq,	double *vloopq,		double *dptavq,
    double *tptavq,	double *betapq,	double *betatq,	double *palphaq,	double *ploplhq);
}


int main(int argc, char *argv[]) {

// ================================================================================
// NEW - Read in setup parameters from setup.data input file
	ifstream fin_disp;  fin_disp.open("InputData/setup.data");
	fin_disp.width(10); fin_disp.precision(2);
	fin_disp.setf(ios_base::uppercase, ios_base::scientific);
	int itxt_disp = 28;  // length of setup.data line
	
	int toggle_coe;      fin_disp.seekg(24*itxt_disp); fin_disp >> toggle_coe;
	int toggle_system;   fin_disp.seekg(27*itxt_disp); fin_disp >> toggle_system;
	int toggle_sic;      fin_disp.seekg(15*itxt_disp); fin_disp >> toggle_sic;
	int toggle_dcll;   	 fin_disp.seekg(17*itxt_disp); fin_disp >> toggle_dcll;
	int toggle_manifold; fin_disp.seekg(31*itxt_disp); fin_disp >> toggle_manifold;

	fin_disp.close();

// ================================================================================
// Define display mode - "COE" or "system"
extern string dispmode; // "coe" is for production mode for scanning ALL input points
dispmode = "coe";
if (toggle_system == 1) { dispmode = "system"; } // "system" is for detailed analysis and geometry on a SINGLE point

// ================================================================================
// Define blanket to use - SCLL or DCLL
extern string blanketmode; string plasmaswitch = "Plasma_AT"; // default to sic
blanketmode = "sic";
if (toggle_dcll == 1) { blanketmode = "dcll"; plasmaswitch = "Plasma_DCLL"; }

// ================================================================================
// Select certain components
extern string manifold;
manifold = "off"; // default is off
if (toggle_manifold == 1) { manifold = "on"; }

cout << " Blanket = " << blanketmode << ", mode = " << dispmode << 
		", manifolds are " << manifold << ", "; // printout

// ================================================================================

// LCC - why are these three commented out below??? 
// Maybe this is from an all-in-one program that goes from inpar files directly thru 
// physics, engineering, and the Matlab engineering filters.  
// Or perhaps it is commented out to allow the full range to run and not constrain any parameter.
// Perhaps they have not been deleted when the program was migrated to C++

/*   double rq[900], aspq[900], btq[900], aipq[900], betanq[900],
     qmhdq[900], qcylq[900], deltaq[900], anq[900], atq[900], fracgwq[900],
     Qq[900], akq[900], hfactorq[900], taujq[900], taueq[900], tpoteq[900],
     xneq[900], tempaveq[900], tflatq[900], pltohq[900], xnwallq[900],
     pbremq[900], fbsq[900], etacdq[900], powcdq[900], pauxq[900], pcyclq[900],
     pohmq[900], plineq[900], betafq[900], zeffq[900], tratq[900], dratq[900],
     powcd1q[900], powcd2q[900], powcd3q[900], powcd4q[900], fimp1q[900],
     fimp2q[900], fimp3q[900], temp0q[900], fcdq[900], fniq[900], fheq[900],
     fdtq[900], wtotq[900], dfcsq[900], vloopq[900], dptavq[900], tptavq[900],
     betapq[900], betatq[900], palphaq[900], ploplhq[900];                      */

// Run the Physics Module
/*   sys_(&rq[0], &aspq[0], &btq[0], &aipq[0], &betanq[0],
     &qmhdq[0], &qcylq[0], &deltaq[0], &anq[0], &atq[0], &fracgwq[0],
     &Qq[0], &akq[0], &hfactorq[0], &taujq[0], &taueq[0], &tpoteq[0],
     &xneq[0], &tempaveq[0], &tflatq[0], &pltohq[0], &xnwallq[0],
     &pbremq[0], &fbsq[0], &etacdq[0], &powcdq[0], &pauxq[0], &pcyclq[0],
     &pohmq[0], &plineq[0], &betafq[0], &zeffq[0], &tratq[0], &dratq[0],
     &powcd1q[0], &powcd2q[0], &powcd3q[0], &powcd4q[0], &fimp1q[0],
     &fimp2q[0], &fimp3q[0], &temp0q[0], &fcdq[0], &fniq[0], &fheq[0],
     &fdtq[0], &wtotq[0], &dfcsq[0], &vloopq[0], &dptavq[0], &tptavq[0],
     &betapq[0], &betatq[0], &palphaq[0], &ploplhq[0]);                      */

// Run the Engineering Filters
/*   sysengr_(&rq[0], &aspq[0], &btq[0], &aipq[0], &betanq[0],
     &qmhdq[0], &qcylq[0], &deltaq[0], &anq[0], &atq[0], &fracgwq[0],
     &Qq[0], &akq[0], &hfactorq[0], &taujq[0], &taueq[0], &tpoteq[0],
     &xneq[0], &tempaveq[0], &tflatq[0], &pltohq[0], &xnwallq[0],
     &pbremq[0], &fbsq[0], &etacdq[0], &powcdq[0], &pauxq[0], &pcyclq[0],
     &pohmq[0], &plineq[0], &betafq[0], &zeffq[0], &tratq[0], &dratq[0],
     &powcd1q[0], &powcd2q[0], &powcd3q[0], &powcd4q[0], &fimp1q[0],
     &fimp2q[0], &fimp3q[0], &temp0q[0], &fcdq[0], &fniq[0], &fheq[0],
     &fdtq[0], &wtotq[0], &dfcsq[0], &vloopq[0], &dptavq[0], &tptavq[0],
     &betapq[0], &betatq[0], &palphaq[0], &ploplhq[0]);                      */


// Read input data from physics module
int len = 1100000;  // LCC changed from 300k to 900k to 1.1M, thought it might be limiting scan
double (*chuck)[100];
chuck = new double [len][100]; // "Chuck" = Chuck Kessel

// Translate the code external input argument into the data point number, which is integer
 int ndp;

 if( strlen(argv[1]) == 1)
 ndp = 1*(argv[1][0] - '0');
 
 if( strlen(argv[1]) == 2)
 ndp = 10*(argv[1][0] - '0') + 1*(argv[1][1] - '0');

 if( strlen(argv[1]) == 3)
 ndp = 100*(argv[1][0] - '0') + 10*(argv[1][1] - '0') + 1*(argv[1][2] - '0');
 
 if( strlen(argv[1]) == 4)
 ndp = 1000*(argv[1][0] - '0') + 100*(argv[1][1] - '0') +
 10*(argv[1][2] - '0') + 1*(argv[1][3] - '0');
 
 if( strlen(argv[1]) == 5)
 ndp = 10000*(argv[1][0] - '0') + 1000*(argv[1][1] - '0') +
 100*(argv[1][2] - '0') + 10*(argv[1][3] - '0') + 1*(argv[1][4] - '0');
 
 if( strlen(argv[1]) == 6)
 ndp = 100000*(argv[1][0] - '0') + 10000*(argv[1][1] - '0') +
 1000*(argv[1][2] - '0') + 100*(argv[1][3] - '0') + 10*(argv[1][4] - '0') +
 1*(argv[1][5] - '0');

// LCC added 1M statement, replicating above, not sure it this helps with files > 1M points...
 if( strlen(argv[1]) == 7)
 ndp = 1000000*(argv[1][0] - '0') + 100000*(argv[1][1] - '0') +
 10000*(argv[1][2] - '0') + 1000*(argv[1][3] - '0') + 100*(argv[1][4] - '0') +
 10*(argv[1][5] - '0') + 1*(argv[1][6] - '0');


// Opens physics module's output file to read the data point it, specified in runaries.sh
	ifstream fin;  
	if (blanketmode == "dcll") { 
		   fin.open("InputData/sysoutdcll.data"); cout << "database = sysoutdcll.data\n"; }
	else { fin.open("InputData/sysoutsic.data");  cout << "database = sysoutsic.data\n"; }
	fin.width(10); fin.precision(2);
	fin.setf(ios_base::uppercase, ios_base::scientific);
	fin.seekg(ndp*661); // length of one row

  for(int j=0; j<55; j++)
  { fin >> chuck[ndp][j]; }
cout << " Point in database: " << ndp << "\n";

system_generator(chuck, ndp, dispmode, blanketmode, plasmaswitch, manifold);
 
} // end of "main" function


//=======================================================================================
// .... .... System Generator Code ... ....
void system_generator(double (*chuck)[100], int ndp, 				string dispmode, 
					  string blanketmode, 	string plasmaswitch, 	string manifold)
{

// Define Parts, links to DP ~46
   double plasma_coeff[7]; 			// added SOL to list of "plasma_params"

// For generic radial build
   double inboard_generic_build[26];
   double outboard_generic_build[28];
   double vertical_generic_build[22];

// Plasma Coefficients from input file
    plasma_coeff[0] = chuck[ndp][0]; 				// Major radius
	plasma_coeff[1] = chuck[ndp][0]/chuck[ndp][1];	// Minor radius
	plasma_coeff[2] = chuck[ndp][12];				// Elongation
	plasma_coeff[3] = chuck[ndp][7];				// Triangularity

// NEW - Read the plasma squareness coefficients from powerflow.data input file
	ifstream fin_plasma;  fin_plasma.open("InputData/powerflow.data");
	fin_plasma.width(10); fin_plasma.precision(2);
	fin_plasma.setf(ios_base::uppercase, ios_base::scientific);
	int itxt_pl = 57;  // length of first line of powerflow.data
	int inum_pl = 13;  // the number called is 13 letters long
	int iskip_pl = itxt_pl + inum_pl; // combination goes to the next parameter
	
	// CK says to use ~0 for now since there are no hard parameters specified for the plasma
	fin_plasma.seekg(47*iskip_pl+itxt_pl+1); fin_plasma >> plasma_coeff[4]; // outer sq +1.0e-6
	fin_plasma.seekg(48*iskip_pl+itxt_pl+1); fin_plasma >> plasma_coeff[5]; // inner sq -1.0e-6

	fin_plasma.close();


// Physics Parameters
	double aip   = chuck[ndp][3];  // plasma current
	double qmhd  = chuck[ndp][5];  // q95, plasma edge safety factor
	double nwl   = chuck[ndp][21]; // average neutron wall load (NWLave) from physics module
	double betap = chuck[ndp][51]; // plasma poloidal beta


//=======================================================================================
// Read Radial Build - SCLL and DCLL, radial build from plasma edge out to PF coil

// See "20090423 - El-Guebaly Blankets.pdf" for radial builds + compositions
// Also LCC "ASC Legend 24.xls" compositions tab
// NEW - generic radial builds implemented, "builds_generic_sic.data"
// Blankets: "dcll" = April 2009, "dcll2" = July 2011, "sic" = AT-type

  ifstream fin_build;
  if (manifold == "on" && blanketmode == "dcll")
  								  fin_build.open("InputData/builds_generic_dcll2.data");
  else if (blanketmode == "dcll") fin_build.open("InputData/builds_generic_dcll.data");
  else 					    	  fin_build.open("InputData/builds_generic_sic.data");
								 
  fin_build.setf(ios_base::uppercase, ios_base::scientific);
  fin_build.width(10);	fin_build.precision(2);
  int line = 28;  // length of builds_generic input file line

// Find size of file
/* long begin, end;
 begin = fin_build.tellg();
 fin_build.seekg (0, ios::end);
 end = fin_build.tellg();
 cout << "generic build file size is: " << (end-begin) << " bytes.\n";
*/

// ... inboard ... (generic starts at line 6)
  if (blanketmode == "dcll") {
	for(int j=0; j<26; j++) { fin_build.seekg((j*2+6)*line); fin_build >> inboard_generic_build[j]; } }
//cout << j << " " << inboard_generic_build[j] << "\n";
  else {
	for(int j=0; j<26; j++) { fin_build.seekg((j*2+6)*line); fin_build >> inboard_generic_build[j]; } }


// ... outboard ... (generic starts at line 62)
  if (blanketmode == "dcll") {
	for(int j=0; j<28; j++) { fin_build.seekg((j*2+62)*line); fin_build >> outboard_generic_build[j]; } }
  else {
  	for(int j=0; j<28; j++) { fin_build.seekg((j*2+62)*line); fin_build >> outboard_generic_build[j]; } }


// ... vertical ... (generic starts at line 122)
  if (blanketmode == "dcll") {
	for(int j=0; j<22; j++) { fin_build.seekg((j*2+122)*line); fin_build >> vertical_generic_build[j]; } }
  else {
  	for(int j=0; j<22; j++) { fin_build.seekg((j*2+122)*line); fin_build >> vertical_generic_build[j]; } }

fin_build.close();



// Modify NWL with plasma surface area to include SOL, copied from DP.ccp ~240
	plasma_coeff[6] = inboard_generic_build[0]; // SOL thickness, done for plasma_coeff above
/*	double SOL = inboard_generic_build[0];
    double P_neutron = 4.0*chuck[ndp][53]; // = 4.0*P_alpha = Pn
    const double  pi = 4.0*atan(1.0);
	double A = chuck[ndp][1];
	double k = chuck[ndp][12];
	double R = chuck[ndp][0];
	double A_plasma_surface = 4.0*pi*pi*R*(R/A + SOL) *
					   		  sqrt((1.0 + pow(k, 2.0))/ 2.0); // using elliptical formula

	//double A_plasma_surface = (*Plasma).surface_area("outer"); // too early to calculate
*/	
	// FOR NOW, use nwl from physics, more accurate than using elliptical formulas
    //double NWL_ave_w_SOL = P_neutron/A_plasma_surface;
	//nwl  = NWL_ave_w_SOL; // reclassify physics nwl to include area with SOL



//=======================================================================================
// SCALE HT shield thickness normalized to NWL = 3.3
// Include the effect of the NWL_ave on the HT shield thickness (reduces it). 
// Only inboard HT shield and top/bottom repl HT shields are affected.
// NOTE: in C++ "log" is actually the natural log (ln base e).  Log base 10 is "log10()"

// HT Shield Scaling Factors
	double ht_sf;
	if (blanketmode == "dcll") ht_sf = 7.5e-2*log(nwl/3.3); // DCLL
	else					   ht_sf = 6.7e-2*log(nwl/3.3); // SCLL

	if  (inboard_generic_build[8] == 0.0 ) // do not modify if thickness is zero
		 inboard_generic_build[8] = inboard_generic_build[8];
	else inboard_generic_build[8] = inboard_generic_build[8] + ht_sf; // use NWL w/ SOL
	
	if  (vertical_generic_build[5] == 0.0 )
		 vertical_generic_build[5] = vertical_generic_build[5];
	else vertical_generic_build[5] = vertical_generic_build[5] + ht_sf; // use NWL w/ SOL
	
// Include the effect of the NWL_ave on the outboard skeleton ring
	if  (outboard_generic_build[15] == 0.0 )
		 outboard_generic_build[15] = outboard_generic_build[15];
	else outboard_generic_build[15] = outboard_generic_build[15] + ht_sf; // use NWL w/ SOL
	
	
//=======================================================================================
// Get Plasma and Radial Builds

// Allows build parameters to be accessed by other functions, e.g. G.cpp
// e.g. (*Machine).inboard_generic_build(1)
DesignPoint Aries_DCLL; DesignPoint Aries_AT; // need to declare before IF statement
if (blanketmode == "dcll") { 
	Aries_DCLL.get_plasma(plasma_coeff);
	Aries_DCLL.get_builds(inboard_generic_build, outboard_generic_build, vertical_generic_build);
	Aries_DCLL.get_physics(nwl, qmhd, aip, betap); }
else { 
	Aries_AT.get_plasma(plasma_coeff);
	Aries_AT.get_builds(inboard_generic_build, outboard_generic_build, vertical_generic_build);
	Aries_AT.get_physics(nwl, qmhd, aip, betap); }

// Initialize the data point filters:
// See below ~600 for addition to data file output
	if (blanketmode == "dcll") { for(int i=0; i<20; i++) Aries_DCLL.get_data_point_filter(0, i); }
	else 					   { for(int i=0; i<20; i++) Aries_AT.get_data_point_filter(0, i); }


//=======================================================================================
// Read in Reference Years for Costing Accounts

// Read reference years from costing.data input file
  ifstream fin_coe;  fin_coe.open("InputData/costing.data");
  fin_coe.width(10); fin_coe.precision(2);
  fin_coe.setf(ios_base::uppercase, ios_base::scientific);
  int itxt_coe = 57;  // length of first line of costing.data
  int inum_coe = 13;  // the number called is 13 letters long
  int iskip_coe = itxt_coe + inum_coe;

  int costBaseYear; fin_coe.seekg(8*iskip_coe + itxt_coe +1); fin_coe >> costBaseYear; // base year reference
  int costOutYear;  fin_coe.seekg(9*iskip_coe + itxt_coe +1); fin_coe >> costOutYear;  // year in use
  // Inflation factor = cost basis year in use/cost reference year, costOutYear/costBaseYear
  fin_coe.close();
 
 
//=======================================================================================
// PRINT OUT inflationfactors.data once (removed print out from Part.cpp ~990)

// Read GDP inflation factors from inflationfactors.data input file
	ifstream fin_gdp;  fin_gdp.open("InputData/inflationfactors.data");
	fin_gdp.width(10); fin_gdp.precision(2);
	fin_gdp.setf(ios_base::uppercase, ios_base::scientific);
	int next = 15;     // next line
	int year[42];      // number of years
	double factor[85]; // number of factors (lines)

// Read the year data from the input file
for(int i=3, j=0; i<86; j++)
    { fin_gdp.seekg(i*next);     fin_gdp >> year[j];   // year
      fin_gdp.seekg((i+1)*next); fin_gdp >> factor[j]; // factor
      i=i+2; 
//cout << year[j] << "	;" << factor[j] << "\n";
}
      
// 12/2010: GDP Current (103010)(2005 basis).xls
// 10/2011: GDP Current (102111)(2005 basis).xls
// 02/2012: GDP Current (012712)(2005 basis).xls
ofstream fileout("./OutputParts/zinflationfactors.data");
fileout << "GDP Inflation Factors used in ASC: \n\n";
fileout << " Factors updated 02/2012 per 'GDP current (012712)(2005 basis).xls' L. Waganer \n";
fileout << " GDP Year basis:	;" << costBaseYear << "\n"; // basis from which table is calculated
fileout << " Year used:	;" << costOutYear << "\n\n";

fileout << "Year	;Inflation Factor \n";
	for(int i=0; i<42; i++)
    	{ fileout 	<< year[i] << "	;" << factor[i] << "\n"; }
    	
fin_gdp.close();


//=======================================================================================
// Read in Material Properties for print out

// Read material properties from materials.data input file
	ifstream fin_mat;  fin_mat.open("InputData/materials.data"); // 1 of 3
	fin_mat.width(10); fin_mat.precision(2);
	fin_mat.setf(ios_base::uppercase, ios_base::scientific);
	int next2 = 60;      // length of line
	string name[22];    // number of materials
	string name2[22];   // number of materials for descriptive name
	string name3[22];   // number of materials for purpose
	double prop[131];   // number of total properties (lines)
	double density[22]; // number of materials

for(int i=0, j=0, c=0; i<22; i++) // 22 materials
    { fin_mat.seekg(j*next2);     fin_mat >> name[i];   // name
      fin_mat.seekg((j+1)*next2); fin_mat >> name2[i];  // descriptive name
      fin_mat.seekg((j+2)*next2); fin_mat >> name3[i];  // purpose
      fin_mat.seekg((j+3)*next2); fin_mat >> prop[c];   // density
	  fin_mat.seekg((j+4)*next2); fin_mat >> prop[c+1]; // cost/kg
	  fin_mat.seekg((j+5)*next2); fin_mat >> prop[c+2]; // cost basis year 
//cout << " name = " << name[i] << "\n";
      c=c+5; j=j+7;
      } // iterate by 5 material properties & 7 lines per material name


// PRINT OUT material properties in one file
// 12/2010: Fabrication Costs per kgRevN (01192011).xls
// 10/2011: Fabrication Costs per kgRevP (10282011).xls
ofstream fileout_b("./OutputParts/zmaterialproperties.data");
fileout_b << "Material Properties used in the Aries Systems Code: \n\n";
fileout_b << "All properties updated 10/2011 'Fabrication Costs per kg RevP (10282011).xls' L. Waganer \n";
fileout_b << "All coil properties provided by L. Bromberg ~2007 'TF, PF Magnet.doc' \n\n";
fileout_b << "Material Name	;Description	;Purpose	;Density (kg/m3)	;$/kg	;Cost basis year \n";
	for(int i=0, c=0; i<22; i++) // 22 materials
    	{ fileout_b	<< name[i] << "	;"
    				<< name2[i] << "	;"
     				<< name3[i] << "	;"
   					<< prop[c] << "	;"
    				<< prop[c+1] << "	;"
    				<< prop[c+2] << "\n";
    				density[i] = prop[c]; // put density in a separate vector to use below
					c=c+5; } // Iterate by 5 material properties
fileout_b << "\n *Helium and liquid metals costs are included in the Special Materials Account 27 since they are added at startup and are not part of the capital costs. \n\n";

fin_mat.close();


//=======================================================================================
//=======================================================================================
// Build each component from the different modules

// All these from "Plasma" to "Upper Vacuum Duct" are from Geometry.cpp file

// 1. Plasma - both (i.e. same design function is called)
Part Plasma_AT;	  
Part Plasma_AT_w_SOL; // needed for separate part class to calc plasma SA with SOL
Part Plasma_DCLL;
if (blanketmode == "dcll") { Plasma_DCLL.get_name("Plasma_DCLL");
	design_plasma_shape(&Aries_DCLL, &Plasma_DCLL);	design_plasma_shape_w_SOL(&Aries_DCLL, &Plasma_AT_w_SOL); 
} // use this name for plasma w/ SOL, even for DCLL
else { Plasma_AT.get_name("Plasma_AT");
	design_plasma_shape(&Aries_AT, &Plasma_AT);
	design_plasma_shape_w_SOL(&Aries_AT, &Plasma_AT_w_SOL); 
}


// 2. Separatrix - both
Part Separatrix;
Separatrix.get_name("Separatrix");
if (blanketmode == "dcll") design_separatrix_shape(&Aries_DCLL, &Separatrix, &Plasma_DCLL);
else                       design_separatrix_shape(&Aries_AT,   &Separatrix, &Plasma_AT);


// 3. First Wall Inboard - both (but DCLL is different - NEED to implement)
Part FirstWallInboard;
FirstWallInboard.get_name("FirstWallInboard");
if (blanketmode == "dcll") design_first_wall_inboard_shape(&Aries_DCLL, &FirstWallInboard, &Plasma_DCLL);
else                       design_first_wall_inboard_shape(&Aries_AT,   &FirstWallInboard, &Plasma_AT);
	cost_first_wall_inboard(&FirstWallInboard, costBaseYear);
	if(dispmode == "system") FirstWallInboard.get_mass();


// 4. First Wall Outboard - both
Part FirstWallOutboard;
FirstWallOutboard.get_name("FirstWallOutboard");
if (blanketmode == "dcll") design_first_wall_outboard_shape(&Aries_DCLL, &FirstWallOutboard, &Plasma_DCLL);
else					   design_first_wall_outboard_shape(&Aries_AT,   &FirstWallOutboard, &Plasma_AT);
	cost_first_wall_outboard(&FirstWallOutboard, costBaseYear);
	if(dispmode == "system") FirstWallOutboard.get_mass();


// 5. Blanket Inboard - both
Part BlanketInboard;
BlanketInboard.get_name("BlanketInboard");
if (blanketmode == "dcll") design_blanket_inboard_shape(&Aries_DCLL, &BlanketInboard, &FirstWallInboard);
else					   design_blanket_inboard_shape(&Aries_AT,   &BlanketInboard, &FirstWallInboard);
	cost_blanket_inboard(&BlanketInboard, costBaseYear);
	if(dispmode == "system") BlanketInboard.get_mass();


// 6. Blanket Outboard - both
Part BlanketOutboard;
BlanketOutboard.get_name("BlanketOutboard");
if (blanketmode == "dcll") design_blanket_outboard_shape(&Aries_DCLL, &BlanketOutboard, &FirstWallOutboard);
else					   design_blanket_outboard_shape(&Aries_AT,   &BlanketOutboard, &FirstWallOutboard);
	cost_blanket_outboard(&BlanketOutboard, costBaseYear);
	if(dispmode == "system") BlanketOutboard.get_mass();


// 7a. Back Wall Inboard - DCLL
Part BackWallInboard; Part BackWallOutboard; Part BlanketII;
if (blanketmode == "dcll") {
BackWallInboard.get_name("BackWallInboard");
design_back_wall_inboard_shape(&Aries_DCLL, &BackWallInboard, &BlanketInboard);
	cost_back_wall_inboard(&BackWallInboard, costBaseYear);
	if(dispmode == "system") BackWallInboard.get_mass();


// 7b. Back Wall Outboard - DCLL
BackWallOutboard.get_name("BackWallOutboard");
design_back_wall_outboard_shape(&Aries_DCLL, &BackWallOutboard, &BlanketOutboard);
	cost_back_wall_outboard(&BackWallOutboard, costBaseYear); 
	if(dispmode == "system") BackWallOutboard.get_mass(); }


// 7c. Blanket II - SCLL
else { BlanketII.get_name("BlanketII");
design_blanket_II_shape(&Aries_AT, &BlanketII, &BlanketOutboard);
	cost_blanket_II(&BlanketII, costBaseYear);
	if(dispmode == "system") BlanketII.get_mass(); }


// 8. Divertor Plates - both
Part Divertor[6];
Divertor[0].get_name("Div1");
Divertor[1].get_name("Div2");
Divertor[2].get_name("Div3");
Divertor[3].get_name("Div4");
Divertor[4].get_name("Div5");
Divertor[5].get_name("Div6");
if (blanketmode == "dcll") design_divertor_shape_dcll(&Aries_DCLL, &Divertor, &Separatrix,
                       	   &BackWallInboard, &FirstWallOutboard, &Plasma_DCLL);
else 					   design_divertor_shape(&Aries_AT, &Divertor, &Separatrix,
                           &BlanketInboard, &BlanketOutboard, &Plasma_AT);
	cost_divertor(&Divertor, costBaseYear);
	if(dispmode == "system") { for(int i=0; i<6; i++) Divertor[i].get_mass(); }


// 9. Replaceable HT Shield - both
Part RepHTShield[2];
RepHTShield[0].get_name("RepHTShieldTop");
RepHTShield[1].get_name("RepHTShieldBottom");
if (blanketmode == "dcll") 
 design_replaceable_ht_shield(&Aries_DCLL, &RepHTShield, &BackWallInboard, &FirstWallOutboard, &Plasma_DCLL); // was &BackWallInboard, &FirstWallOutboard, &Plasma_DCLL
else
 design_replaceable_ht_shield(&Aries_AT, &RepHTShield, &BlanketInboard, &BlanketOutboard, &Plasma_AT);
	cost_replaceable_ht_shield(&RepHTShield, costBaseYear);
	if(dispmode == "system") { for(int i=0; i<2; i++) RepHTShield[i].get_mass(); }


// 10a+b. HT Shield: DCLL maintains it as one unit, SCLL breaks it into IB+OB
Part HTShieldDCLL; Part HTShield[2];
if (blanketmode == "dcll") { HTShieldDCLL.get_name("HTShield");
//if (manifold == "on") ;// do nothing, no HT shield with DCLL2 (gives a bus error...)
     design_ht_shield_dcll(&Aries_DCLL, &HTShieldDCLL, &BackWallInboard, &RepHTShield, &BlanketOutboard, &BackWallOutboard); 
	 cost_ht_shield_dcll(&HTShieldDCLL, costBaseYear);
	if(dispmode == "system") HTShieldDCLL.get_mass(); }
else { HTShield[0].get_name("HTShieldInboard"); 
	   HTShield[1].get_name("HTShieldOutboard");
	design_ht_shield_sic(&Aries_AT, &HTShield, &BlanketInboard, &RepHTShield, &BlanketOutboard, &BlanketII);
	cost_ht_shield_sic(&HTShield, costBaseYear); }
	if(dispmode == "system") { for(int i=0; i<2; i++) HTShield[i].get_mass(); }
	

// 10.5 Skeleton Ring - DCLL and DCLL2
Part SkeletonRing;
if (blanketmode == "dcll") {
	 SkeletonRing.get_name("SkeletonRing");
if (manifold == "on")
	 design_skeleton_ring2(&Aries_DCLL, &SkeletonRing, &BackWallInboard, &RepHTShield, &BlanketOutboard, &BackWallOutboard);
else design_skeleton_ring( &Aries_DCLL, &SkeletonRing, &HTShieldDCLL);
	 cost_skeleton_ring(&SkeletonRing, costBaseYear);
	if(dispmode == "system") SkeletonRing.get_mass(); }
 
 
// 10.6 Manifold - DCLL OB and vert
//if (manifold == "on") {
Part Manifold;
if (blanketmode == "dcll") {
	Manifold.get_name("Manifold");
	design_manifold(&Aries_DCLL, &Manifold, &SkeletonRing);
	cost_manifold(&Manifold, costBaseYear);
 	if(dispmode == "system") Manifold.get_mass(); } 


// 10.7 Manifolding - SCLL
	Part Manifolding[2];
	Manifolding[0].get_name("ManifoldingTop");
	Manifolding[1].get_name("ManifoldingBottom");
if   (blanketmode == "dcll") ; // do nothing yet for DCLL
else design_manifolding(&Aries_AT, &Manifolding, &RepHTShield, &FirstWallOutboard, &BlanketOutboard, &BlanketII);
	cost_manifolding(&Manifolding, costBaseYear);
 	if(dispmode == "system") { for(int i=0; i<2; i++) Manifolding[i].get_mass(); }


// 11a. VacuumVessel - DCLL
Part VacuumVesselDCLL[4];
Part VacuumVessel[4]; // DCLL, SCLL (last 2/4 are just placeholders for SCLL)
if (blanketmode == "dcll") { 
	VacuumVesselDCLL[0].get_name("VacuumVesselInboard");
	VacuumVesselDCLL[1].get_name("VacuumVesselOutboard");
	VacuumVesselDCLL[2].get_name("VacuumVesselTop");
	VacuumVesselDCLL[3].get_name("VacuumVesselBottom");
if (manifold == "on") { 
design_vacuum_vessel_inboard(	&Aries_DCLL, &VacuumVesselDCLL, &Manifold);	      // 1/4
design_vacuum_vessel_outboard(	&Aries_DCLL, &VacuumVesselDCLL, &Manifold);       // 2/4
design_vacuum_vessel_top(		&Aries_DCLL, &VacuumVesselDCLL, &Manifold); }	  // 3/4
else {
design_vacuum_vessel_inboard(	&Aries_DCLL, &VacuumVesselDCLL, &SkeletonRing);	  // 1/4
design_vacuum_vessel_outboard(	&Aries_DCLL, &VacuumVesselDCLL, &SkeletonRing);   // 2/4
design_vacuum_vessel_top(		&Aries_DCLL, &VacuumVesselDCLL, &SkeletonRing); } // 3/4

copy_vacuum_vessel_top_into_bottom_by_symmetry(&VacuumVesselDCLL);            	  // 4/4
	cost_vacuum_vessel_dcll(&VacuumVesselDCLL, costBaseYear);
	if(dispmode == "system") { for(int i=0; i<4; i++) VacuumVesselDCLL[i].get_mass(); } }


// 11b. VacuumVessel - SCLL
else { VacuumVessel[0].get_name("VacuumVesselMain");
	   VacuumVessel[1].get_name("VacuumVesselOutboard");
	design_vacuum_vessel(&Aries_AT, &VacuumVessel, &HTShield, &RepHTShield);
cost_vacuum_vessel(&VacuumVessel, costBaseYear);
if(dispmode == "system") {for(int i=0; i<2; i++) VacuumVessel[i].get_mass(); } }


// 12. TF Coil - both, do we NEED low-temp for SCLL, high-temp for DCLL
// This ARIES-AT coil design is chosen from 7 versions of r-theta magnet algorithms in Geometry.cpp
// Composed of r-theta AND r-z planes
Part TFCoil;
TFCoil.get_name("TFMagnet");
double BT = chuck[ndp][2];

if (blanketmode == "dcll") { 
	design_tf_coil_rtheta_NextStepConservative(&Aries_DCLL, &TFCoil, &Plasma_DCLL, &VacuumVesselDCLL, BT);   
	design_tf_coil_rz_plane(&Aries_DCLL, &TFCoil, &VacuumVesselDCLL);
	cost_tf_coil_conservative(&TFCoil, &Aries_DCLL, costBaseYear); }
else {
	design_tf_coil_rtheta_NextStepConservative(&Aries_AT, &TFCoil, &Plasma_AT, &VacuumVessel, BT); // G.cpp ~1830   
	design_tf_coil_rz_plane(&Aries_AT, &TFCoil, &VacuumVessel);    // G.cpp ~2790
	cost_tf_coil_conservative(&TFCoil, &Aries_AT, costBaseYear); } // Cost.cpp ~147 
if(dispmode == "system" ) { TFCoil.get_number_of_material_components(5); TFCoil.get_mass(); }


// 13. Bucking Cylinder - both
Part BuckingCylinder;
BuckingCylinder.get_name("BuckingCylinder");
if (blanketmode == "dcll") design_bucking_cylinder(&Aries_DCLL, &BuckingCylinder, &TFCoil, &VacuumVesselDCLL, BT);
else					   design_bucking_cylinder(&Aries_AT,   &BuckingCylinder, &TFCoil, &VacuumVessel,     BT);
	cost_bucking_cylinder(&BuckingCylinder, costBaseYear);
	if(dispmode == "system" ) BuckingCylinder.get_mass();


// 14. Upper Cap of TF Coil - both
Part UpperCap;
UpperCap.get_name("UpperCap");
if (blanketmode == "dcll")  design_upper_cap_shape(&Aries_DCLL, &UpperCap, &TFCoil, &VacuumVesselDCLL, BT);  
else						design_upper_cap_shape(&Aries_AT,   &UpperCap, &TFCoil, &VacuumVessel,     BT);
	cost_upper_cap(&UpperCap, costBaseYear);
	if(dispmode == "system" ) UpperCap.get_mass();


 // 15. Lower Cap of TF Coil - both
Part LowerCap;
LowerCap.get_name("LowerCap");
copy_upper_cap_into_lower_cap_shape_by_symmetry(&LowerCap, &UpperCap);   
//design_lower_cap_shape(&Aries_AT, &LowerCap, &TFCoil, &VacuumVessel, BT); // not working
	cost_lower_cap(&LowerCap, costBaseYear);
	if(dispmode == "system" ) LowerCap.get_mass();


// 16. Maintenance Port Window Opening - was SCLL only, now include DCLL
Part MPWindow;
//if (blanketmode == "sic") {
MPWindow.get_name("MaintenancePortOpening");
design_mpw(&Aries_AT, &MPWindow, &TFCoil, &VacuumVessel); // cost is "TF Coil" bc that is what it is a window between
	cost_mpw(&MPWindow, costBaseYear);


// 17. PF Coil
Part PoloFCoil[36];
for(int i=0; i<36; i++)
 { std::ostringstream stm; stm << i; string name = "PFCoil_" + stm.str(); PoloFCoil[i].get_name(name); }
if (blanketmode == "dcll") design_pf_coil(&Aries_DCLL, &PoloFCoil, &BuckingCylinder, &UpperCap, &LowerCap);
else					   design_pf_coil(&Aries_AT,   &PoloFCoil, &BuckingCylinder, &UpperCap, &LowerCap);
	cost_pf_coil(&PoloFCoil, costBaseYear);
	if(dispmode == "system") {for(int i=0; i<36; i++) PoloFCoil[i].get_mass(); }


// 18. Cryogenic Dome or Cryostat
Part CryoDome;
CryoDome.get_name("CryoDome");
if (blanketmode == "dcll") design_cryogenic_dome_shape(&Aries_DCLL, &CryoDome, &TFCoil);
else					   design_cryogenic_dome_shape(&Aries_AT,   &CryoDome, &TFCoil);
	cost_cryogenic_dome(&CryoDome, costBaseYear);
	if(dispmode == "system") CryoDome.get_mass();


// 19. Maintenance Port
Part MaintenancePort;
MaintenancePort.get_name("MaintenancePort");
if (blanketmode == "dcll") design_maintenance_port(&Aries_DCLL, &MaintenancePort, &VacuumVesselDCLL, &CryoDome);
else					   design_maintenance_port(&Aries_AT,   &MaintenancePort, &VacuumVessel, &CryoDome);
	cost_maintenance_port(&MaintenancePort, costBaseYear);


// 20. Upper Vacuum Duct, composed of bottom, top and side walls
Part UpperVVDuct[3];
UpperVVDuct[0].get_name("UpperVVDuctBottomWall");
UpperVVDuct[1].get_name("UpperVVDuctTopWall");
UpperVVDuct[2].get_name("UpperVVDuctSideWall");
if (blanketmode == "dcll") {
	design_upper_vv_duct_bottom_wall(&Aries_DCLL, &UpperVVDuct, &VacuumVesselDCLL, &TFCoil, &CryoDome);
	design_upper_vv_duct_top_wall(   &Aries_DCLL, &UpperVVDuct, &VacuumVesselDCLL, &TFCoil, &CryoDome);
	design_upper_vv_duct_side_wall(  &Aries_DCLL, &UpperVVDuct); }
else {
	design_upper_vv_duct_bottom_wall(&Aries_AT, &UpperVVDuct, &VacuumVessel, &TFCoil, &CryoDome);
	design_upper_vv_duct_top_wall(   &Aries_AT, &UpperVVDuct, &VacuumVessel, &TFCoil, &CryoDome);
	design_upper_vv_duct_side_wall(  &Aries_AT, &UpperVVDuct); }
cost_upper_vv_duct(&UpperVVDuct, costBaseYear);

// All above are called from Geometry.cpp ..................................


//=======================================================================================
// Power Flow Scheme, called from DesignPoint.cpp

// Need only one powerflow scheme, want dcll and scll combined so only uses one DP function.
// Need to add DCLL components to SCLL powerflow then copy Aries_AT to Aries_DCLL - done.
if (blanketmode == "dcll")
Aries_DCLL.get_power_flow(
			&FirstWallInboard,  	&FirstWallOutboard, 	&Divertor, 
			&Plasma_AT, 			chuck, 					ndp, 			
			&BlanketInboard, 		&BlanketOutboard, 		&BlanketII,
			&RepHTShield, 			&HTShield,				&VacuumVessel,		
			&TFCoil, 				&PoloFCoil,				&Plasma_AT_w_SOL,
			&Plasma_DCLL,			&BackWallInboard,		&BackWallOutboard,
			&SkeletonRing,			&VacuumVesselDCLL,		&HTShieldDCLL,
			&Manifold,				&Manifolding);
			// copied Aries_AT.get+power_flow........
else 
Aries_AT.get_power_flow(
			&FirstWallInboard,  	&FirstWallOutboard, 	&Divertor, 
			&Plasma_AT, 			chuck, 					ndp, 			
			&BlanketInboard, 		&BlanketOutboard, 		&BlanketII,
			&RepHTShield, 			&HTShield,				&VacuumVessel,		
			&TFCoil, 				&PoloFCoil,				&Plasma_AT_w_SOL,
			&Plasma_DCLL,			&BackWallInboard,		&BackWallOutboard,
			&SkeletonRing,			&VacuumVesselDCLL,		&HTShieldDCLL,
			&Manifold,				&Manifolding);
			// added for DCLL: Plasma_DCLL on

// ====================================================================================
// Costing

if (blanketmode == "dcll") { 
	Aries_DCLL.get_reactor_building_volume(&CryoDome);
	Aries_DCLL.get_support_structure_volume();
	Aries_DCLL.get_cd_cost(); }
else {
	Aries_AT.get_reactor_building_volume(&CryoDome); // from DP ~54
	Aries_AT.get_support_structure_volume(); // Calc in DP ~63; support structure vol ~20% FPC vol 
	}										 // old - ref for 1.7221e8? was in cm3 (why?), so = 172.21 m3...

// CD parameters are called here from DP.cpp ~1788
Aries_AT.get_cd_cost(); // no manual inputs now, calls from powerflow.data input file
CostingAccount ATCosts; // same
ATCosts.get_cost_base_years(costBaseYear, costOutYear); // transfers cost base years to CA.cpp ~33
ATCosts.get_level_of_safety_assurance(4);  // Set as bogus LSA=4 for cost credit factor input



// DCLL Costing accounts
// SCLL and DCLL accounts should mimic each other and contain ALL the parts that the account needs
// for both designs.  Only difference is first term &Aries_DCLL or &Aries_AT.
if (blanketmode == "dcll") { 
ATCosts.get_account_20(		&Aries_DCLL);
ATCosts.get_account_21(		&Aries_DCLL);
ATCosts.get_account_22_1(  	&Aries_DCLL, 		&FirstWallInboard, 	&FirstWallOutboard,
                           	&BlanketInboard,	&BlanketOutboard, 	&BlanketII,
                           	&RepHTShield, 		&HTShield, 			&TFCoil, 
                           	&MPWindow,      	&PoloFCoil, 		&BuckingCylinder, 
                           	&UpperCap, 			&LowerCap,          &VacuumVessel, 
                           	&Divertor,			&BackWallInboard,	&BackWallOutboard,
                           	&HTShieldDCLL,		&SkeletonRing,		&VacuumVesselDCLL,
                           	&Manifold);
ATCosts.get_account_22_2_1(	&Aries_DCLL,		&TFCoil, 			&BuckingCylinder, 
						   	&UpperCap,       	&LowerCap, 			&CryoDome);
ATCosts.get_account_22_2(	&Aries_DCLL, 		&TFCoil, 			&BuckingCylinder, 
							&UpperCap,		 	&LowerCap, 			&PoloFCoil, 		
							&CryoDome);ATCosts.get_account_22_3_1(	&Aries_DCLL);
ATCosts.get_account_22_3_2(	&Aries_DCLL);
ATCosts.get_account_22_3_3(	&Aries_DCLL);
ATCosts.get_account_22_3_4(	&Aries_DCLL);
ATCosts.get_account_22_3(	&Aries_DCLL);						   	
ATCosts.get_account_22_4(	&Aries_DCLL, 		&VacuumVessel,	 	&MaintenancePort, 
							&UpperVVDuct,		&VacuumVesselDCLL); // second &VacuumVesselDCLL is placeholder for DCLL
ATCosts.get_account_22_5(	&Aries_DCLL);
ATCosts.get_account_22_6(	&Aries_DCLL);
ATCosts.get_account_22_7(	&Aries_DCLL);
ATCosts.get_account_22_8(	&Aries_DCLL);
ATCosts.get_account_22_9(	&Aries_DCLL);
ATCosts.get_account_22_10(	&Aries_DCLL);
ATCosts.get_account_22_11(	&Aries_DCLL);
ATCosts.get_account_22(		&Aries_DCLL); // note, main account must be after sub accounts
ATCosts.get_account_23(		&Aries_DCLL);
ATCosts.get_account_24(		&Aries_DCLL);
ATCosts.get_account_25(		&Aries_DCLL);
ATCosts.get_account_26(		&Aries_DCLL);
ATCosts.get_account_27(		&Aries_DCLL, 		&FirstWallInboard, 	&FirstWallOutboard,
                       		&BlanketInboard, 	&BlanketOutboard, 	&BlanketII,
                       		&Divertor, 			&RepHTShield, 		&HTShield,
                       		&BackWallInboard,	&BackWallOutboard,	&HTShieldDCLL,
                       		&SkeletonRing,		&Manifold);                       		
} // end DCLL


// SCLL Costing accounts
else {
ATCosts.get_account_20(	   	&Aries_AT);
ATCosts.get_account_21(    	&Aries_AT);
ATCosts.get_account_22_1(  	&Aries_AT, 			&FirstWallInboard, 	&FirstWallOutboard,
                          	&BlanketInboard, 	&BlanketOutboard, 	&BlanketII,
                          	&RepHTShield, 		&HTShield, 			&TFCoil, 
                           	&MPWindow,       	&PoloFCoil, 		&BuckingCylinder, 
                           	&UpperCap, 			&LowerCap,          &VacuumVessel, 
                           	&Divertor,			&BackWallInboard,	&BackWallOutboard,
                           	&HTShieldDCLL,		&SkeletonRing,		&VacuumVesselDCLL,
                           	&Manifold); 
                           // added for DCLL: BackWallInboard on
ATCosts.get_account_22_2_1( &Aries_AT, 			&TFCoil, 			&BuckingCylinder, 
						    &UpperCap,       	&LowerCap, 			&CryoDome);
ATCosts.get_account_22_2(   &Aries_AT, 			&TFCoil, 			&BuckingCylinder, 
						    &UpperCap,       	&LowerCap, 			&PoloFCoil, 
						    &CryoDome);
ATCosts.get_account_22_3_1(	&Aries_AT);
ATCosts.get_account_22_3_2(	&Aries_AT);
ATCosts.get_account_22_3_3(	&Aries_AT);
ATCosts.get_account_22_3_4(	&Aries_AT);
ATCosts.get_account_22_3(	&Aries_AT);						   
ATCosts.get_account_22_4(	&Aries_AT, 			&VacuumVessel, 		&MaintenancePort, 	
							&UpperVVDuct,		&VacuumVesselDCLL); // second &VacuumVessel is placeholder for DCLL
ATCosts.get_account_22_5(	&Aries_AT); 
ATCosts.get_account_22_6(	&Aries_AT);
ATCosts.get_account_22_7(	&Aries_AT);
ATCosts.get_account_22_8(	&Aries_AT);
ATCosts.get_account_22_9(	&Aries_AT);
ATCosts.get_account_22_10(	&Aries_AT);
ATCosts.get_account_22_11(	&Aries_AT);
ATCosts.get_account_22(	   	&Aries_AT);
ATCosts.get_account_23(		&Aries_AT);
ATCosts.get_account_24(		&Aries_AT);
ATCosts.get_account_25(		&Aries_AT);
ATCosts.get_account_26(		&Aries_AT);
ATCosts.get_account_27(		&Aries_AT, 			&FirstWallInboard, 	&FirstWallOutboard,
                       		&BlanketInboard, 	&BlanketOutboard, 	&BlanketII,
                       		&Divertor, 			&RepHTShield, 		&HTShield,
                       		&BackWallInboard,	&BackWallOutboard,	&HTShieldDCLL,
                       		&SkeletonRing,		&Manifold); 
} // end SCLL

// SCLL and DCLL shared accounts                    
ATCosts.get_account_90();
ATCosts.get_account_91();
ATCosts.get_account_92();
ATCosts.get_account_93();
ATCosts.get_account_94();
ATCosts.get_account_95();
ATCosts.get_account_96();
ATCosts.get_account_97();
ATCosts.get_account_98();
ATCosts.get_account_99();


// Cost of Electicity, from CostingAccount.cpp
	if (blanketmode == "dcll") ATCosts.get_cost_of_electricity(&Aries_DCLL);
	else 					   ATCosts.get_cost_of_electricity(&Aries_AT);


// ================================================================================
// Fill the input/output vector with new data:

if (blanketmode == "dcll") {
   chuck[ndp][21] = Aries_DCLL.neutron_wall_load("average");
   chuck[ndp][55] = Plasma_DCLL.surface_area("outer");
   chuck[ndp][56] = Aries_DCLL.tf_magnet_max_field();
   chuck[ndp][57] = Aries_DCLL.maximum_heat_flux_on_divertor("inboard");
   chuck[ndp][58] = Aries_DCLL.maximum_heat_flux_on_divertor("outboard");
   chuck[ndp][59] = Aries_DCLL.fw_heat_flux("maximum");
   chuck[ndp][60] = Aries_DCLL.net_electrical_power();
   chuck[ndp][61] = Aries_DCLL.gross_electrical_power() - Aries_DCLL.net_electrical_power();
   chuck[ndp][62] = Aries_DCLL.thermal_power();
   chuck[ndp][63] = ATCosts.cost_of_electricity();
   chuck[ndp][64] = ndp; }

else {
   chuck[ndp][21] = Aries_AT.neutron_wall_load("average"); // modifies [21] from MW/m2 to W/m2...what implication does this have??
   chuck[ndp][55] = Plasma_AT.surface_area("outer");
   chuck[ndp][56] = Aries_AT.tf_magnet_max_field();
   chuck[ndp][57] = Aries_AT.maximum_heat_flux_on_divertor("inboard");
   chuck[ndp][58] = Aries_AT.maximum_heat_flux_on_divertor("outboard");
   chuck[ndp][59] = Aries_AT.fw_heat_flux("maximum");
   chuck[ndp][60] = Aries_AT.net_electrical_power();
   chuck[ndp][61] = Aries_AT.gross_electrical_power() - Aries_AT.net_electrical_power();
   chuck[ndp][62] = Aries_AT.thermal_power();
   chuck[ndp][63] = ATCosts.cost_of_electricity(); 
   chuck[ndp][64] = ndp; } // data point number


// Add the data point filters to the output file from [65] to [83].
// If the filter fails, then it seems as if the failed filter just puts its number in its filter 
// slot.  E.g. if Pnelec fails (DP ~677), then a "7" is inserted into column 64+7=71. This could 
// have had a filtering application prior, but it seems better to simply filter the value as we 
// see fit with an external program (Matlab) rather than hard-wiring filter values in the code.
for(int i=1; i<20; i++) 
	if (blanketmode == "dcll") chuck[ndp][64+i] = Aries_DCLL.data_point_filter(i);
	else					   chuck[ndp][64+i] = Aries_AT.data_point_filter(i);


// ================================================================================
// "SYSTEM" mode
// Display detailed output system data for a single data point

// need - problem is here with dcll, bus error after this, make separate Aries_DCLL.show???

if(dispmode == "system") { 
	double OutputArray[100]; 
  	for(int j=0; j<84; j++) // number of parameters in the sysout file
  	{ OutputArray[j] = chuck[ndp][j]; } // puts the sysout from physics file into new output array

// Shows the detailed printout info from DP.cpp
// Show sas components for both SCLL and DCLL. Aries_AT and Aries_DCLL are each able to go
// into "DesignPoint::show()" function in DP ~1040
// Add DCLL components to SCLL show then copy Aries_AT to Aries_DCLL 
 
if (blanketmode == "dcll")	   
 Aries_DCLL.show( &Plasma_AT, 		&Separatrix, 		&FirstWallInboard,
               &FirstWallOutboard, 	&BlanketInboard,	&BlanketOutboard, 
               &BlanketII, 			&Divertor, 			&RepHTShield,
               &HTShield, 			&VacuumVessel, 		&TFCoil, 
               &MPWindow, 			&BuckingCylinder,   &UpperCap, 
               &LowerCap,           &PoloFCoil, 		&CryoDome,
               &MaintenancePort, 	&UpperVVDuct, 		&Plasma_AT_w_SOL,
               &OutputArray,		&Plasma_DCLL,		&SkeletonRing,
               &BackWallInboard,	&BackWallOutboard,	&VacuumVesselDCLL,
               &HTShieldDCLL,		&Manifold,			&Manifolding);
               // copy of Aries_AT.show
else 
Aries_AT.show( &Plasma_AT, 			&Separatrix, 		&FirstWallInboard,
               &FirstWallOutboard, 	&BlanketInboard,	&BlanketOutboard, 
               &BlanketII, 			&Divertor, 			&RepHTShield,
               &HTShield, 			&VacuumVessel, 		&TFCoil, 
               &MPWindow, 			&BuckingCylinder,   &UpperCap, 
               &LowerCap,           &PoloFCoil, 		&CryoDome,
               &MaintenancePort, 	&UpperVVDuct, 		&Plasma_AT_w_SOL,
               &OutputArray,		&Plasma_DCLL,		&SkeletonRing,
               &BackWallInboard,	&BackWallOutboard,	&VacuumVesselDCLL,
               &HTShieldDCLL,		&Manifold,			&Manifolding);		
               // added for DCLL: Plasma_DCLL - on
// NOTE the special blanketmode case uses &,* unlike power_flow. Also note that * is for Part
 
 ATCosts.show(); // shows costs
 } // end dispmode system
 
 
// ================================================================================
// "COE" mode
// Fill the output file \COE Maps\SysoutFinalSiC.data when scanning all data points

if(dispmode == "coe" && blanketmode == "sic") {
	ofstream fileout_a("./COEMaps/SysoutFinalSiC.data",  ios_base::out | ios_base::app);  
  for(int j=0; j<84; j++)
  { fileout_a.width(14);
    fileout_a.precision(6);
    fileout_a.setf(ios_base::uppercase, ios_base::scientific);
    fileout_a << chuck[ndp][j]; }
    fileout_a << "\n";
    fileout_a.close(); 
   } // end dispmode coe

if(dispmode == "coe" && blanketmode == "dcll") {
	ofstream fileout_a("./COEMaps/SysoutFinalDCLL.data", ios_base::out | ios_base::app);  
  for(int j=0; j<84; j++)
  { fileout_a.width(14);
    fileout_a.precision(6);
    fileout_a.setf(ios_base::uppercase, ios_base::scientific);
    fileout_a << chuck[ndp][j]; }
    fileout_a << "\n";
    fileout_a.close(); 
   } // end dispmode coe


delete [] chuck;
} // end "system_generator" function