// Member Functions of the Class "DesignPoint"
// SiC blanket
// LCC cleaned up formatting 9/13/10

#include <cstdlib>
#include <iostream>
#include <iomanip>  // LCC added to manipulate output file format, such as using setw()
#include <fstream>
#include <string>
#include <math.h>
#include <time.h> // LCC for date/time printout on the top of System.data

using std::string;
using namespace std;

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

// Declare global variables from Aries.cpp
extern string blanketmode, dispmode, manifold;

// ====================================================================================
// Get information from other modules, see also line 1537

void DesignPoint::get_parts(int np)
{ cout << "Got parts. \n"; }

void DesignPoint::get_plasma(double params[7]) // from Aries.cpp ~290
{ for(int i=0; i<7; i++) plasma_params[i] = params[i]; }


void DesignPoint::get_xpoint(double x, double y)
{ x_point_radius = x;  x_point_height = y; }

void DesignPoint::get_tf_coil_radii(double r1, double r2)
{ tf_coil_inboard_radius  = r1;  tf_coil_outboard_radius = r2; }

void DesignPoint::get_maximum_magnetic_field_at_tf_coil(double field)
{ tf_coil_max_magnetic_field = field; }

void DesignPoint::get_number_of_tf_cable_turns(int turns) 
{ number_of_tf_cable_turns = turns; }


// ====================================================================================
// Get radial builds

// links to line "get.build" in Aries.cpp
void DesignPoint::get_builds(double inbg[26], double outbg[28], double vertg[22])
{ for(int i=0; i<28; i++) // max mem location is OB side (27)
  	{ rad_gen_inb[i]  = inbg[i];
  	  rad_gen_outb[i] = outbg[i];
  	  top_gen_bld[i]  = vertg[i]; } } // eof


// ====================================================================================
// Reactor Building Volume

// Uses exterior building dimensions, called by CA.cpp ~123, Accnt 21.2
void DesignPoint::get_reactor_building_volume(Part *Core) // Core = CryoDome
{   double ksihi  = (*Core).bounding_box(2, "outer"); // OB core width (x)
    double etahi  = (*Core).bounding_box(3, "outer"); // vert core height (y)
    double cATfit = 138698.0/181655.0; 	// calibration factor to fit the ARIES-AT = 0.76352426, KEEP .0!
    // Fit Prometheus-HI vol = 181655 (ref?), AT vol = 138698 m3 (pub. lit.)
    vol_reactor_bldg = (pow((ksihi + 9.0), 2.0)*6.0*etahi + 1.55e5)*cATfit;  // reference?
//cout << "vol_reactor_bldg (DP.cpp) = " << vol_reactor_bldg << "\n"; // 141,358 m3
//cout << "ksihi = " << ksihi << "\n"; 
//cout << "etahi = " << etahi << "\n";
} // eof


// ====================================================================================
// Support Structure Volume

// called by Aries.cpp ~527, was a constant ~172 m3, should be 20% of FPC vol
void DesignPoint::get_support_structure_volume()
{	double FPC_vol = fpc_vol(); // call DP.cpp ~780
	vol_str = 0.20*FPC_vol;
//cout << "Support structure volume = " << vol_str << "\n";
//cout << "FPC volume2 = " << FPC_vol << "\n";
}
	
	
// ====================================================================================
// ====================================================================================
// Power Flow Module

void DesignPoint::get_power_flow(
		Part * FirstWallInb, 		Part * FirstWallOutb,		Part (*Div)[6],
		Part * Plasma,				double (*inputdata)[100],	int npt, 			
		Part * BlanketInboard,		Part * BlanketOutboard,		Part * BlanketII,
		Part (*RepHTShield)[2],		Part (*HTShield)[2],		Part (*VacuumVessel)[4],	
		Part * TFCoil,				Part (*PoloFCoil)[36],		Part * Plasma_AT_w_SOL,
		Part * Plasma_DCLL,			Part * BackWallInboard,		Part * BackWallOutboard,
		Part * SkeletonRing,		Part (*VacuumVesselDCLL)[4],Part * HTShieldDCLL,
		Part * Manifold,			Part (*Manifolding)[2]) 
		// note: inputdata is called "chuck", added for DCLL: Plasma_DCLL on
{
// Open the input file "Powerflow.data" to read the power flow parameters:
// LCC "inputdata" is the sysout.data input file array
  ifstream fin_pow;  fin_pow.open("InputData/powerflow.data");
  fin_pow.width(10); fin_pow.precision(2);
  fin_pow.setf(ios_base::uppercase, ios_base::scientific);
  int itxt_pow = 57; 	// LCC length of line describing the parameter
  int inum_pow = 13; 	// LCC length of containing the value
  int iskip_pow = itxt_pow + inum_pow;
  const double pi = 4.0*atan(1.0);


// Power Conversion Factors..........................
// 0. fraction of loss of alpha particles = 0 for tokamak
	double f_alpha_loss; fin_pow.seekg(0*iskip_pow + itxt_pow +1); fin_pow >> f_alpha_loss; 	

// 1. fraction of alpha loss to divertor, 50%
	double f_alpha_div; fin_pow.seekg(1*iskip_pow + itxt_pow +1); fin_pow >> f_alpha_div; 	

// 2. fraction of edge radiation from particles (original value was 0.75)
	double f_div_rad; fin_pow.seekg(2*iskip_pow + itxt_pow +1); fin_pow >> f_div_rad;		

// 3. fraction of edge radiation to divertor
	double f_rad_edge; fin_pow.seekg(3*iskip_pow + itxt_pow +1); fin_pow >> f_rad_edge;

// Fraction of edge radiation to FWBS = (1.0 - f_rad_edge)
// IB/OB split of edge radiation to FWBS is dependent on FW area:
	double f_rad_edge_blanket_ib = (*FirstWallInb).surface_area("outer")/
								  ((*FirstWallInb).surface_area("outer")+(*FirstWallOutb).surface_area("outer"));
	double f_rad_edge_blanket_ob = (1.0 - f_rad_edge_blanket_ib);
//cout << f_rad_edge_blanket_ib << "\n";

// Neutron Energy Multiplication (all M = 1.1 for now) ..................
// 4. FW/blanket IB
	double f_nmul_fw_inboard; fin_pow.seekg(4*iskip_pow + itxt_pow +1); fin_pow >> f_nmul_fw_inboard; 

// 5. FW/blanket OB
	double f_nmul_fw_outboard; fin_pow.seekg(5*iskip_pow + itxt_pow +1); fin_pow >> f_nmul_fw_outboard; 	

// 6. shield IB
	double f_nmul_shield_inboard; fin_pow.seekg(6*iskip_pow + itxt_pow +1); fin_pow >> f_nmul_shield_inboard; 		

// 7. shield OB
	double f_nmul_shield_outboard; fin_pow.seekg(7*iskip_pow + itxt_pow +1); fin_pow >> f_nmul_shield_outboard; 
	
// 8. divertor plates
	double f_nmul_divertor; fin_pow.seekg(8*iskip_pow + itxt_pow +1); fin_pow >> f_nmul_divertor; 		

// 9. blanket behind divertor
	double f_nmul_div_blanket; fin_pow.seekg(9*iskip_pow + itxt_pow +1); fin_pow >> f_nmul_div_blanket; 	

// 10. shield behind divertor
	double f_nmul_div_shield; fin_pow.seekg(10*iskip_pow + itxt_pow +1); fin_pow >> f_nmul_div_shield; 	

// 57-59. LT components behind divertor, behind FWBS IB & OB
//	double f_nmul_div_LT;   fin_pow.seekg(57*iskip_pow + itxt_pow +1); fin_pow >> f_nmul_div_LT; 	
//	double f_nmul_fw_IB_LT; fin_pow.seekg(58*iskip_pow + itxt_pow +1); fin_pow >> f_nmul_fw_IB_LT; 	
//	double f_nmul_fw_OB_LT; fin_pow.seekg(59*iskip_pow + itxt_pow +1); fin_pow >> f_nmul_fw_OB_LT; 	


// Efficiency of pumps
// 11, 49. Efficiency of the He pump/blower
	double eta_pump_he, eta_pump_he_recovered;
	fin_pow.seekg(11*iskip_pow + itxt_pow +1); fin_pow >> eta_pump_he; // 0.9
	fin_pow.seekg(49*iskip_pow + itxt_pow +1); fin_pow >> eta_pump_he_recovered; // 0.9

// 34. Efficiency of the Liquid Metal MHD pump 
// (not very practical - Siegfried says will not get much head or P > 2 Mpa)
	double eta_pump_lm_mhd; // 0.4
	double eta_pump_lm_mhd_recovered = 0.8; // Suggested by M. Tillack, N. Morley
    fin_pow.seekg(34*iskip_pow + itxt_pow +1); fin_pow >> eta_pump_lm_mhd;

// 35, 50. Efficiency of the Liquid Metal mechanical pump
	double eta_pump_lm_mech, eta_pump_lm_mech_recovered; 
    fin_pow.seekg(35*iskip_pow + itxt_pow +1); fin_pow >> eta_pump_lm_mech; // 0.9
	fin_pow.seekg(50*iskip_pow + itxt_pow +1); fin_pow >> eta_pump_lm_mech_recovered; // 0.9

	fin_pow.close();
// ...............................................................................
// LCC is this the max number of "fin_pow.seekg" lines read???  Stops reading here
// REINSTATE read input file lines with "fin_po" from "fin_pow"
  ifstream fin_po;  fin_po.open("InputData/powerflow.data");
  fin_po.width(10); fin_po.precision(2);
  fin_po.setf(ios_base::uppercase, ios_base::scientific);
  int itxt_po = 57; int inum_po = 13; int iskip_po = itxt_po + inum_po;    


// ====================================================================================
// Fusion Power
// Power from neutron particles
    double P_alpha   = inputdata[npt][53]*1.0e6; // alpha power contained in plasma, Watts
    double P_neutron = 4.0*P_alpha;
    double P_fusion  = P_alpha + P_neutron; // P_fusion = 5*P_alpha

// 12. Get generic CD efficiency from powerflow.data input file
    double eta_cd_generic; fin_po.seekg(12*iskip_po + itxt_po +1); fin_po >> eta_cd_generic;
    double P_cd[4], eta_cd[4];
    int nb =   0; // sets up an array so P_cd[0]=nb, P_cd[1]=lh, ...
    int lh =   1;
    int icrf = 2;
    int ec =   3;

// Powers are taken from the physics module
// These are the defining parameters for the CD source

    P_cd[nb] =   0.0*inputdata[npt][34]*1.0e6; // These parameters from physics file are not in use
    P_cd[lh] =   0.0*inputdata[npt][35]*1.0e6;
    P_cd[icrf] = 0.0*inputdata[npt][36]*1.0e6;
    P_cd[ec] =   0.0*inputdata[npt][37]*1.0e6;

// Efficiencies are taken from published literature (from powerflow.data input file)
// 13-16. Additional heating sources include NB, LH, ICRF, EC
    fin_po.seekg(13*iskip_po + itxt_po +1); fin_po >> eta_cd[nb];
    fin_po.seekg(14*iskip_po + itxt_po +1); fin_po >> eta_cd[lh];
    fin_po.seekg(15*iskip_po + itxt_po +1); fin_po >> eta_cd[icrf];
    fin_po.seekg(16*iskip_po + itxt_po +1); fin_po >> eta_cd[ec];

// Current Drive Power
    double P_cd_generic = inputdata[npt][25]*1.0e6;
    double P_cd_loss 	= (1.0 - eta_cd_generic)*P_cd_generic/eta_cd_generic;
      
         
// 17. Power to Auxiliary Functions:
// was set at at 50 MW in input file, 2 MW cryo separate, need detailed list of what is included
// L. El-Guebaly recommends 65 MW = 30 misc + 35 MW cryo load (Dec 2010), scaling from ITER
    double P_aux_func; fin_po.seekg(17*iskip_po + itxt_po +1); fin_po >> P_aux_func;

// 18. Cryogenic Magnet Power:
// "use 2 MW for low temperature superconductor and 0.5 MW for high temperature SC"
    double P_cryo; fin_po.seekg(18*iskip_po + itxt_po +1); fin_po >> P_cryo; 

// 19. Plasma Heating Power,  [26] = aux. power into plasma for power balance
    double eta_plasma_heating; fin_po.seekg(19*iskip_po + itxt_po +1); fin_po >> eta_plasma_heating;
	double P_aux_for_power_balance = (inputdata[npt][26])*1.0e6; // includes CD power
    double P_plasma_heating 	   = P_aux_for_power_balance - P_cd_generic; // in MW
	double P_plasma_heating_loss   = (1.0 - eta_plasma_heating)*P_plasma_heating/eta_plasma_heating;

// Power Radiated from Chamber:
// Bremsstralhlung radiation + Cyclotron power + Line radiation
    double P_rad_chamb = (inputdata[npt][22] + inputdata[npt][27] + inputdata[npt][29])*1.0e6;


// ====================================================================================
// Plasma Surface Area

// OLD way was not right, SA should NOT include divertor. 
// C. Kessel: BEST to use plasma + SOL surface area = 4*pi^2*R*(a+SOL)*sqrt((1+k^2)/2)
// CHANGED "A_plasma_facing" to "A_plasma_surface"
// = ~497 m2 using perim method, 512 using eqn, 549 using eqn w/ SOL 10 cm (see Misc tab on spreadsheet)

/* Ellipical equations - good for estimation, but use CK's more accurate discretized formulas in Geometry ~30
	double SOL = rad_gen_inb[0];     // SOL thickness from input file, 10cm
	double A   = inputdata[npt][1];  // Aspect ratio from physics module output, A=4.0
	double k   = inputdata[npt][12]; // kappa, plasma elongation 1.8-2.2
	double R   = inputdata[npt][0];  // plasma major radius or =((*Plasma).bounding_box(2,"outer") + (*Plasma).bounding_box(0,"outer"))/2.0;
	double plasma_major_R = R;
*/

// Perimeter method
	//double A_plasma_surface = (*Plasma).surface_area("plasma_perimeter")*2.0*pi*plasma_major_R; 

// Elliptical method
	//double A_plasma_surface = 4.0*pi*pi*R * (R / A + SOL) * sqrt((1.0 + pow(k, 2.0))/ 2.0); 

// Discrete Functions Method: from Geometry ~100 "Plasma Geometry WITH SOL"
double A_plasma_surface, A_plasma_surface_w_SOL, A_first_wall;
if (blanketmode == "dcll")
	A_plasma_surface 	   	= (*Plasma_DCLL).surface_area("outer"); // does NOT include SOL
else {
	A_plasma_surface 	   	= (*Plasma).surface_area("outer"); // does NOT include SOL
	A_plasma_surface_w_SOL 	= (*Plasma_AT_w_SOL).surface_area("outer"); } // includes SOL

    A_first_wall 			= ((*FirstWallInb).surface_area("outer") +
                               (*FirstWallOutb).surface_area("outer"))/2.0*16.0;					   

//cout << "A_plasma_surface = " << A_plasma_surface << "\n";
//cout << "Area first wall = " << A_first_wall << "\n";


// ====================================================================================
// Geometry factor for divertor, proportions power to FW+blanket/divertor based on area fractions
// ... 16x segmentation factor, then divertors are multiplied x2 to cover upper and lower

/*OLD double f_geo_div = 16*2*((*Div)[0].surface_area("outer") +
                              (*Div)[1].surface_area("outer") +
                              (*Div)[2].surface_area("outer"))/A_plasma_facing;  // 0.274144

// M.Tillack, better estimate is theta/pi, waiting for input from L. El-Guebaly
	double f_geo_div = 0.10942;  // Hand-calc from strawmen CAD as theta into divertor / pi */


// 44-46. Neutron power distribution from plasma given from L. El-Guebaly, read from powerflow.data
// "For A=4, 10% of power is deposited in divertor, 25% in IB blanket, and 65% in OB.  
    double f_n_div;        fin_po.seekg(44*iskip_po + itxt_po +1); fin_po >> f_n_div;
    double f_n_blanket_ib; fin_po.seekg(45*iskip_po + itxt_po +1); fin_po >> f_n_blanket_ib;
    double f_n_blanket_ob; fin_po.seekg(46*iskip_po + itxt_po +1); fin_po >> f_n_blanket_ob;

// NEW: 51-53. Radiated power distribution from plasma, read from powerflow.data
    double f_rad_div;        fin_po.seekg(51*iskip_po + itxt_po +1); fin_po >> f_rad_div;
    double f_rad_blanket_ib; fin_po.seekg(52*iskip_po + itxt_po +1); fin_po >> f_rad_blanket_ib;
    double f_rad_blanket_ob; fin_po.seekg(53*iskip_po + itxt_po +1); fin_po >> f_rad_blanket_ob;


// ====================================================================================
// Average Neutron Wall Flux (aka "average neutron wall load")

// C.K. - neutron wall load = 4/5*Pfusion/plasma_SA = P_neutron/A_plasma_surface
    double flux_neutron_wall_Ave 	   = P_neutron/A_plasma_surface; // A_plasma_surface calc above ~240
    double flux_neutron_wall_Ave_w_SOL = P_neutron/A_plasma_surface_w_SOL;    
//cout << "NWL ave = " << flux_neutron_wall_Ave << "\n";
//cout << "P_neutron = " << P_neutron << "\n";


// Maximum Neutron Wall Load Table using Peaking Factor; linear interpolation inbetween
// Laila updated SiC 4/4/11, at A=4, PF from 1.5 to 1.68, 5/27/11 to 1.55
   double tab_A[2], tab_PF[2];
if(blanketmode == "dcll"){
      tab_A [0] = 1.6,  tab_A [1] = 4.00;
      tab_PF[0] = 1.94, tab_PF[1] = 1.50; } // at A=4, was 1.71
else{ tab_A [0] = 1.7,  tab_A [1] = 4.00;
      tab_PF[0] = 2.0,  tab_PF[1] = 1.55; }

    double plasma_aspect_ratio = inputdata[npt][1];
    double peaking_factor =  tab_PF[0] +
           					(tab_PF[1] - tab_PF[0])*
           					(plasma_aspect_ratio - tab_A[0])/
           					(tab_A[1]-tab_A[0]);
           					
    // Laila wants limit < 5 MW/m2 but no need - we can filter out later...
    double flux_neutron_wall_Max = flux_neutron_wall_Ave*peaking_factor; 
    get_neutron_wall_loads(flux_neutron_wall_Ave, flux_neutron_wall_Max);
//cout << "NWL max = " << flux_neutron_wall_Max << "\n";


// Power from alpha particles:
// LCC changed from P_plasma to P_alpha since loss is only from P_alpha
// If there is any alpha loss, it would go into the FW and be recovered as heat
    double P_alpha_loss = f_alpha_loss * P_alpha; // = 0 for a tokamak
	double P_transport  = P_alpha - P_alpha_loss;
// LCC changed P_alpha to P_transport
    double P_plasma     = P_transport + P_aux_for_power_balance; // in MW
    double P_particle   = P_plasma - P_alpha_loss - P_rad_chamb;       
//cout << "P_alpha_loss 		= " << P_alpha_loss/1.0e6 << "\n";
//cout << "P_transport 		= " << P_transport/1.0e6 << "\n";
//cout << "P_plasma 		= " << P_plasma/1.0e6 << "\n";
//cout << "P_particle 		= " << P_particle/1.0e6 << "\n";

// Edge radiation power:
    double P_rad_edge = f_div_rad * P_particle;      

// Edge conducted power:
    double P_cond_div = (1.0 - f_div_rad) * P_particle;      


// Heat Flux to First Wall IB + OB (Average and Maximum): 
// LCC changed (1.0-f_geo_div) to f_rad_blanket_ib+f_rad_blanket_ob
// This has increased P_rad_first_wall as well as flux_heat_wall to ~0.298 MW
    double P_rad_first_wall = (f_rad_blanket_ib+f_rad_blanket_ob)	*P_rad_chamb + 
    						  (1.0 - f_rad_edge)					*P_rad_edge;

// Average heat flux on the wall:
    double  flux_heat_wall_Ave = P_rad_first_wall/A_first_wall;

// Maximum heat flux on the wall:
    double flux_heat_wall_Max = 1.25*flux_heat_wall_Ave;
    get_fw_heat_fluxes(flux_heat_wall_Ave, flux_heat_wall_Max);
    

// ====================================================================================
// ====================================================================================
// BLANKET PUMPING POWERS for SiC and DCLL blankets

// LCC removed old methods which divided pumping power 50:50 between blanket:divertor
// on 5/9/11 and placed in "removedcode.cpp"                         

// see "Tillack - SiC Blanket Pumping Power.xls" 3/7/2011
// NOTES: IB FW/B and PbLi ONLY; OB will be much lower due to lower field strength and more space
// USE ONLY P_rad_first_wall_IB, NWL_ave_IB, P_blanket_IB.

// Note: must declare all variables before entering the IF statement
double P_rad_first_wall_IB, A_FW_IB, flux_neutron_wall_IB, P_blanket_sic, 
	   P_blanket, P_blanket_dcll, P_pump_blanket, P_pump_blanket_loss;

// SiC Blanket PUMPING POWER Curve
if(blanketmode == "sic") {
    P_rad_first_wall_IB  = (f_rad_blanket_ib) * P_rad_chamb + 
    					   (1.0 - f_rad_edge) * P_rad_edge  * f_rad_edge_blanket_ib;
    					   
	A_FW_IB 			 = 16.0/2.0*( (*FirstWallInb).surface_area("outer") ); // ~133 m3
	flux_neutron_wall_IB = P_neutron*f_n_blanket_ib/A_FW_IB; // ~2.88 MW/m2. IB+OBave ~2.79
	
	P_blanket_sic 		 = (1.0 - f_alpha_div) * P_alpha_loss + // = 0
						   f_n_blanket_ib     * P_neutron * f_nmul_fw_inboard + 
					       P_rad_first_wall_IB; // = ~460 MW IB ONLY (1172 MW OB)

	P_pump_blanket = (0.43840*(flux_neutron_wall_Ave/1.0e6) + 0.00036) *
					 P_blanket_sic/100.0; // ~5.6 MW
	P_blanket = P_blanket_sic;
	P_pump_blanket_loss = (1.0 - eta_pump_lm_mech)*P_pump_blanket/eta_pump_lm_mech;

// Stay within the range of data given in spreadsheet (filter range is for both SiC and DCLL)
	if(flux_neutron_wall_Ave > 4.0e6  || flux_neutron_wall_Ave < 2.0e6)
	{ cout << "NWL (" << flux_neutron_wall_Ave/1.0e6 << ") is outside the tabluated range 2-4 MW/m2, P_pump_blanket unreliable! \n";
	  P_pump_blanket = 1.0e30; }
	  
//cout << "P_rad_first_wall_IB 	= "	<< P_rad_first_wall_IB/1.0e6 << "\n";
//cout << "P_blanket_sic 		= "	<< P_blanket_sic/1.0e6 << "\n";
}

// ***********OR**************OR*****************OR***************OR*******************
// DCLL Blanket PUMPING POWER Curve

// see "Tillack - DCLL Blanket Pumping Power.xls" 1/14/2011
// INCLUDES: He AND Pb-Li MHD pumping powers are both combined, as well as IB+OB, need to separate...
// blanket thermal power = neutron power*neutron energy multipl + first wall surface heating

//	not using: P_blanket = flux_neutron_wall_Ave *f_nmul_fw_outboard *A_first_wall + 
//			    		   flux_heat_wall_Ave    *A_first_wall; // = ~ 1.8 GW

if(blanketmode == "dcll"){
// Leave area out of P_blanket. Per Laila, 0.9 is P_neutron into FWBS (f_n_blanket_ib + f_n_blanket_ob)
	P_blanket_dcll = (1.0 - f_alpha_div) * P_alpha_loss // = 0
				   + (f_n_blanket_ib) 	 * P_neutron * f_nmul_fw_inboard
				   + (f_n_blanket_ob) 	 * P_neutron * f_nmul_fw_outboard
				   + P_rad_first_wall; // = 1.63 GW, both IB+OB

	P_pump_blanket = (0.03602*pow(flux_neutron_wall_Ave/1.0e6,2) +
					  0.03056  * (flux_neutron_wall_Ave/1.0e6) + 
					  0.38359) *  P_blanket_dcll/100.0; // ~12 MW
					  
// Need coolant fractional power split for DCLL blanket pump losses (LM pump vs He-blower)
// Mark's DCLL pumping power spreadsheet has PbLi 62.3%, He 37.7%
	double f_pump_PbLi = 0.623;
	double f_pump_He   = 0.377;
	P_pump_blanket_loss = ((1.0 - eta_pump_lm_mech)*f_pump_PbLi
						 + (1.0 - eta_pump_he)	   *f_pump_He) 
						 * P_pump_blanket/((eta_pump_he+eta_pump_lm_mech)/2);

// Stay within the range of data given in spreadsheet (filter range is for both SiC and DCLL)
	if(flux_neutron_wall_Ave > 4.0e6  || flux_neutron_wall_Ave < 2.0e6)
	{ cout << "NWL (" << flux_neutron_wall_Ave/1.0e6 << ") is outside the tabluated range 2-4 MW/m2, P_pump_blanket unreliable! \n";
	  P_pump_blanket = 1.0e30; }

//cout << "P_blanket_dcll w/o pp	= "	<< P_blanket_dcll/1.0e6 << "\n";
}

//cout << "NWL max 			= " << flux_neutron_wall_Max << "\n"; 
//cout << "Max heat flux FW 	= " << flux_heat_wall_Max << "\n"; 
//cout << "P_rad_first_wall 	= "	<< P_rad_first_wall/1.0e6 << "\n";
//cout << "P_pump_blanket 		= "	<< P_pump_blanket/1.0e6 << "\n";


// ====================================================================================
// ====================================================================================
// Maximum Heat Flux on Divertor

// from Chuck Kessel's "sysengr.f" systems code
// LCC Max heat flux on divertor calculation moved here from line 550:

/* Definitions
c..........power to divertor
c...sym  		=  1.0 for up-down symmetric, 0.0 for up-down asymmetric
c...facdiv  	=  factor for splitting power if up-down symmetric or not
c...adivo 	 	=  area in outboard divertor slot for absorbing radiated power
c...adivi  		=  area in inboard divertor slot for absorbing radiated power
c...fdivrad  	=  total radiated power fraction in divertor
c...fdivcondo  	=  fraction of conducted power to divertor going to the outboard
c...fdivcondi  	=  fraction of conducted power to divertor going to the inboard
c...fdivrado  	=  fraction of radiated power to divertor going to the outboard
c...fdivradi  	=  fraction of radiated power to divertor going to the inboard
c...fluxexpo  	=  poloidal magnetic flux expansion on outboard from midplane to divertor target
c...fluxexpi  	=  poloidal magnetic flux expansion on inboard from midplane to divertor target   */


// 27. Checks symmetry factor "switch" from powerflow.data, 1 for symmetric, 0 otherwise
   double sym; fin_po.seekg(27*iskip_po + itxt_po +1); fin_po >> sym; // use "po", not "pow"
   
   double facdiv; // divertor symmetry
   if(sym == 1.0) { facdiv = 0.65; } // use DCLL = 0.60?
   else           { facdiv = 1.00; }
 
   double Rmaj = plasma_major_radius();
   double Asp  = plasma_major_radius()/plasma_minor_radius();
      
   double adivo = 2.*pi*(Rmaj -((Rmaj/Asp)/2.0))*2.0*(Rmaj/(Asp*2.0));  // edited w/ CK 1/2011
   double adivi = 2.*pi*(Rmaj-(Rmaj/Asp))*2.0*(Rmaj/(Asp*4.0));  		// edited w/ CK 1/2011
   
// 28-33. Fraction of conducted/radiated power to div, IB, OB
   double fdivcondo; fin_po.seekg(28*iskip_po + itxt_po +1); fin_po >> fdivcondo; // 80%
   double fdivcondi; fin_po.seekg(29*iskip_po + itxt_po +1); fin_po >> fdivcondi; // 20%
   double fdivrado;  fin_po.seekg(30*iskip_po + itxt_po +1); fin_po >> fdivrado;  // 80%
   double fdivradi;  fin_po.seekg(31*iskip_po + itxt_po +1); fin_po >> fdivradi;  // 20%
   double fluxexpo;  fin_po.seekg(32*iskip_po + itxt_po +1); fin_po >> fluxexpo;  // 10
   double fluxexpi;  fin_po.seekg(33*iskip_po + itxt_po +1); fin_po >> fluxexpi;  // 10

/* Definitions
c...pdiv        =  total power going to divertors
c...pdivrad     =  total power radiated in divertors
c...pdiv2       =  total power conducted to divertors
c...pdiv2o      =  power conducted to outboard divertor
c...pdiv2i      =  power conducted to inboard divertor
c...pdivrado    =  power radiated in outboard divertor
c...pdivradi    =  power radiated in inboard divertor
c...pdivtoto    =  conducted + radiated power in outboard divertor
c...pdivtoti    =  conducted + radiated power in inboard divertor
c...den         =  density factor in power scrape-off formula
c...xlam        =  power scrape-off layer width at plasma midplane
c...qpeakcondo  =  peak conducted power heat flux to outboard divertor
c...qpeakcondi  =  peak conducted power heat flux to inboard divertor
c...qpeakrado   =  peak radiated power heat flux to outboard divertor
c...qpeakradi   =  peak radiated power heat flux to inboard divertor  */

   double P_div 		= (P_plasma - P_rad_first_wall)/1.e6;  // the following are all in MW!!
   double P_div_rad 	= f_div_rad*P_div;
   double P_div_2 		= P_div - P_div_rad;
   double P_div_2_o 	= facdiv*fdivcondo*P_div_2;
   double P_div_2_i 	= facdiv*fdivcondi*P_div_2;
   double P_div_rad_o 	= facdiv*fdivrado*P_div_rad;
   double P_div_rad_i 	= facdiv*fdivradi*P_div_rad;

   double f_greenwald 	= inputdata[npt][10]; // greenwald fraction
   double Bt			= inputdata[npt][2];  // toroidal magnetic field at plasma major radius
   double Aip 			= inputdata[npt][3];  // plasma current
   double Qmhd 			= inputdata[npt][5];  // qmhd

   double den  = f_greenwald*(Aip/pi/pow(Rmaj/Asp, 2.0));

   double xlam = 2.9e-2*2.5*pow(Qmhd, 0.75)*pow(den, 0.15)/(pow(P_div, 0.4)*Bt); // power scrape off width
//  7/16/12 modification per F.Najmabadi to use only outboard power (80% of total)
//   double xlam = 2.9e-2*2.5*pow(Qmhd, 0.75)*pow(den, 0.15)/(pow(0.8*P_div, 0.4)*Bt);
// removed after long debate with Kessel 1/10/13

   double Q_peak_cond_o = P_div_2_o/(2.0*pi*(Rmaj-((Rmaj/Asp)/2.0))*fluxexpo*xlam);  // all in MW
   double Q_peak_cond_i = P_div_2_i/(2.0*pi*(Rmaj-(Rmaj/Asp))      *fluxexpi*xlam);
   double Q_peak_rad_o  = P_div_rad_o/adivo;
   double Q_peak_rad_i  = P_div_rad_i/adivi;
   double Q_peak_tot_o  = Q_peak_cond_o + Q_peak_rad_o;
   double Q_peak_tot_i  = Q_peak_cond_i + Q_peak_rad_i;

// Get's peak divertor loads from a simple comparison
   get_divertor_loads(Q_peak_tot_i*1.e6, Q_peak_tot_o*1.e6);
   double flux_divertor_Max = max(Q_peak_tot_o, Q_peak_tot_i)*1.e6;   // this one in Watts

//   cout << "xlam =" << xlam << "\n";
//   cout << "den =" << den << "\n";
//   cout << "P_div =" << P_div << "\n";
//   cout << "P_plasma =" << P_plasma << "\n";
//   cout << "P_rad_first_wall =" << P_rad_first_wall << "\n";
//   cout << "P_rad_chamb =" << P_rad_chamb << "\n";
//   cout << "P_rad_edge =" << P_rad_edge << "\n";
//   cout << "P_div_rad_o =" << P_div_rad_o << "\n";
//   cout << "Q_peak_cond_o =" << Q_peak_cond_o << "\n";
//   cout << "Q_peak_rad_o =" << Q_peak_rad_o << "\n";
//   cout << "Q_peak_tot_o =" << Q_peak_tot_o << "\n";
//   cout << "Max Divertor heat flux = "  << flux_divertor_Max << "\n";

// Filter max heat flux at divertor > 8 MW/m2
 if(flux_divertor_Max > 8.0e6) { get_data_point_filter(6, 6); }


// ====================================================================================
// ====================================================================================
// Divertor Pumping Power: calls line ~1164 for pumping_power_divertor() function

// Note, use P_thermal_DIVERTOR, NOT P_thermal TOTAL

	double P_thermal_divertor = P_alpha_loss *f_alpha_div // = 0
							  + P_neutron	 *f_n_div *f_nmul_divertor
							  + P_rad_edge	 *f_rad_edge
							  + P_rad_chamb	 *f_rad_div
							  + P_cond_div;               // in Watts!
// This section was added by M. Tillack on 7/16/12 to account for tailoring
// of the divertor to the local heat flux.  Two regions are considered:
// a peak flux region (p) and an off-peak flux region (op).
// Both inboard and outboard have these regions, so there are 4 regions.

// We assume a value for the fraction of the total power in the peak vs. off-peak regions.
// (Once we have actual profiles, we can do this right).  Assume 50% for now:
    double fracdivpeak = 0.5;

// we also need to distribute the total P_div 80/20 to the outboard/inboard
    double Pdiv_inboard = 0.2*P_div;              // inboard divertor power, MW
    double Pdiv_outboard = 0.8*P_div;             // outboard divertor power, MW

// we know the peak heat flux (e.g. Q_peak_tot_o), but not the off-peak heat flux.  
// we need to compute it.

// first find the area over which the peak and off-peak are applied
    double adivo_p = Pdiv_outboard * fracdivpeak / Q_peak_tot_o;
    double adivo_op = adivo - adivo_p;
    double adivi_p = Pdiv_inboard * fracdivpeak / Q_peak_tot_i;
    double adivi_op = adivi - adivi_p;

// the off-peak heat flux times the off-peak area equals the off-peak power
    double Q_offpeak_o = Pdiv_outboard*(1.0-fracdivpeak) / adivo_op;
    double Q_offpeak_i = Pdiv_inboard*(1.0-fracdivpeak) / adivi_op;

// next we call the "pumping_power_divertor" subroutines four times,
// to get the pumping power fraction in the four regions, and add them up
    double P_pump_divertor = (
       pumping_power_divertor(Q_peak_tot_i*1.e6) * 0.2 * fracdivpeak       // inboard peak
     + pumping_power_divertor(Q_peak_tot_o*1.e6) * 0.8 * fracdivpeak       // outboard peak
     + pumping_power_divertor(Q_offpeak_i*1.e6) * 0.2 * (1.0-fracdivpeak)  // inboard off-peak
     + pumping_power_divertor(Q_offpeak_o*1.e6) * 0.8 * (1.0-fracdivpeak)  // outboard off-peak
     ) * P_thermal_divertor / 100.0;

/*          
// old code, when we assumed the same heat flux everywhere for pumping power
// P_pump_divertor is the total pumping power in all divertor plates
	double P_pump_divertor = pumping_power_divertor(flux_divertor_Max)
							 *P_thermal_divertor/100.0; // (subroutine uses percents)
*/

	double P_pump_divertor_loss = (1.0 - eta_pump_he)*P_pump_divertor/eta_pump_he;
	// P_pump_divertor/eta_pump_he is the wall plug power req.

// Stay within the range of data given in spreadsheet
	if(flux_divertor_Max > 15.0e6  || flux_divertor_Max < 5.0e6)
	{ cout << "flux_divertor_Max (" << flux_divertor_Max/1.0e6 << ") is outside the tabluated range 5-15 MW/mw, P_pump_divertor unreliable! \n";
	  P_pump_divertor = 1.0e30; }

//cout << "flux_divertor_Max 	= " << flux_divertor_Max/1.0e6 << "\n";
//cout << "P_thermal_divertor w/o pp = " << P_thermal_divertor/1.0e6 << "\n"; // pp not incl
//cout << "P_pump_divertor 	= " << P_pump_divertor/1.0e6 << "\n"; 


// ====================================================================================
// ====================================================================================
// Power absorbed by First Wall / Blanket / Shield: (FWBS)

    double P_FWBSGroup[2][2][2]; // array of powers deposited in FWBS group    
    int neutron    = 0; // i's
    int condRad    = 1;
    int fw_blanket = 0; // j's
    int shield     = 1;
	int outboard   = 0; // k's
    int inboard    = 1;
    
// <o> Incoming Powers:
// i.e. this is what the FWBS sees and is "absorbed"
// LCC changed (1.0-f_geo_div) to f_n_blanket_ib+f_n_blanket_ob (= 90%)

//	double P_neutron_fwbs_ib	  = f_n_blanket_ib*P_neutron;
//	double P_neutron_fwbs_ob	  = f_n_blanket_ob*P_neutron;

// Find LM mechanical pump's thermal contribution to the FWBS Group
// Essentially = wall plug power*eta_pump = P_pump_blanket
	double P_lm_pump_recovered 	  = P_pump_blanket;

    double P_neutron_FWBSGroup_in = (f_n_blanket_ib + f_n_blanket_ob)	* P_neutron;    
    double P_neutron_FWBSGroup_ib = f_n_blanket_ib						* P_neutron;    
    double P_neutron_FWBSGroup_ob = f_n_blanket_ob						* P_neutron;    
    
    double P_condRad_FWBSGroup_in = (1.0 - f_alpha_div)					* P_alpha_loss // = 0
                                  + (f_rad_blanket_ib+f_rad_blanket_ob)	* P_rad_chamb
                                  + (1.0 - f_rad_edge) 					* P_rad_edge;

//cout << " P_neutron_FWBSGroup_ib = " << P_neutron_FWBSGroup_ib/1.0e6 << "\n"; 
//cout << " P_neutron_FWBSGroup_ob = " << P_neutron_FWBSGroup_ob/1.0e6 << "\n"; 
//cout << " P_condRad_FWBSGroup_ib = " << P_condRad_FWBSGroup_ib/1.0e6 << "\n"; 
//cout << " P_condRad_FWBSGroup_ob = " << P_condRad_FWBSGroup_ob/1.0e6 << "\n"; 
    

// 20, 21. Split the incoming/absorbed conducted and radiated power into OB/IB fractions
// These values are recommended by Laila El-Guebaly, based on neutron fractions, 60%, 40%
// Are these for P_rad_chamb? These are superceeded by PFi.45-46, 25%, 65% * P_neutron
//    double f_ncf_outboard; fin_po.seekg(20*iskip_po + itxt_po +1); fin_po >> f_ncf_outboard; 	
//    double f_ncf_inboard;  fin_po.seekg(21*iskip_po + itxt_po +1); fin_po >> f_ncf_inboard;    	


// Conduction + radiation powers are split here. They do not get deposited in the shield,
// only in the FW and blanket, which are the surface layers exposed to plasma.
// Redundant! LCC changed f_ncf_outboard, inboard to f_rad_blanket_ib, ob
    P_FWBSGroup[condRad][fw_blanket][outboard] = f_rad_blanket_ob    * P_rad_chamb
                                  			   + (1.0 - f_rad_edge)  * P_rad_edge
								  			   * f_rad_edge_blanket_ob;
    P_FWBSGroup[condRad][fw_blanket][inboard]  = (1.0 - f_alpha_div) * P_alpha_loss // = 0, some fraction needs to go to ib/ob
                                  			   + f_rad_blanket_ib	 * P_rad_chamb
                                  			   + (1.0 - f_rad_edge)  * P_rad_edge
                                  			   * f_rad_edge_blanket_ib;
    P_FWBSGroup[condRad][shield][outboard]     = 0.0;
    P_FWBSGroup[condRad][shield][inboard]      = 0.0;


// 22-23, 55-56. Split the IB & OB powers from neutrons into FW/blanket and
// shield fractions. Include the effect of neutron energy multiplication (1.1). 90%, 10%
	double f_fw_blanket; fin_po.seekg(22*iskip_po + itxt_po +1); fin_po >> f_fw_blanket;
	double f_shield;     fin_po.seekg(23*iskip_po + itxt_po +1); fin_po >> f_shield;
//	double f_LT_IB;	     fin_po.seekg(55*iskip_po + itxt_po +1); fin_po >> f_LT_IB;
//	double f_LT_OB;	     fin_po.seekg(56*iskip_po + itxt_po +1); fin_po >> f_LT_OB;


 P_FWBSGroup[neutron][fw_blanket][outboard] = f_fw_blanket*				// 0.9
                                              f_nmul_fw_outboard*		// 1.1
                                              P_neutron_FWBSGroup_ob;
                                                 
 P_FWBSGroup[neutron][fw_blanket][inboard] =  f_fw_blanket*				// 0.9
                                              f_nmul_fw_inboard*		// 1.1
                                              P_neutron_FWBSGroup_ib;

 P_FWBSGroup[neutron][shield][outboard] =     f_shield*					// 0.1
                                              f_nmul_shield_outboard*	// 1.1
                                              P_neutron_FWBSGroup_ob;

 P_FWBSGroup[neutron][shield][inboard] =      f_shield*					// 0.1
                                              f_nmul_shield_inboard*	// 1.1
                                              P_neutron_FWBSGroup_ib;
// Add LT VV (~150-200C) to this P_FWBSGroup IB, OB??

/*cout << " P_FWBSGroup[neutron][fw_blanket][outboard] 	= " << P_FWBSGroup[neutron][fw_blanket][outboard]/1.0e6 << "\n"; 
cout << " P_FWBSGroup[neutron][fw_blanket][inboard] 	= " << P_FWBSGroup[neutron][fw_blanket][inboard]/1.0e6 << "\n"; 
cout << " P_FWBSGroup[neutron][shield][outboard] 	= " << P_FWBSGroup[neutron][shield][outboard]/1.0e6 << "\n"; 
cout << " P_FWBSGroup[neutron][shield][inboard] 		= " << P_FWBSGroup[neutron][shield][inboard]/1.0e6 << "\n"; 
cout << " P_FWBSGroup[condRad][fw_blanket][outboard] 	= " << P_FWBSGroup[condRad][fw_blanket][outboard]/1.0e6 << "\n"; 
cout << " P_FWBSGroup[condRad][fw_blanket][inboard] 	= " << P_FWBSGroup[condRad][fw_blanket][inboard]/1.0e6 << "\n"; 
cout << " P_FWBSGroup[condRad][shield][outboard] 	= " << P_FWBSGroup[condRad][shield][outboard]/1.0e6 << "\n"; 
cout << " P_FWBSGroup[condRad][shield][inboard] 		= " << P_FWBSGroup[condRad][shield][inboard]/1.0e6 << "\n"; 
*/

// Note: talk to Laila if there is a need to add the low temperature components.
// Less than 1% of neutron power goes to these components, according to her. (ZD)
// To be implemented when Laila can better tell fraction to LT components.

// FED paper 2006 - Nuclear Performance Assessment of ARIES-AT "The small
// heat leakage to theVV(<1%) is dumped as a low-grade heat."

// Low-temp components (e.g. vacuum vessel)
// Need M?
//	double LT_fwbs_loss = (f_LT_IB*P_neutron_fwbs_ib*f_nmul_fw_IB_LT)
//						+ (f_LT_OB*P_neutron_fwbs_ob*f_nmul_fw_OB_LT);


// ====================================================================================
// ====================================================================================
// Power absorbed by the divertor:

    double P_DivGroup[3][4];  // array of powers deposited in Divertor group
    int divPlates  = 0; // j's
    int divBlanket = 1;
    int divShield  = 2;
	// i's are already defined above in fwbs as neutron and condRad

// <o> Incoming Powers:
// LCC changed f_geo_div to f_n_div

// Find He pump's thermal contribution to the Divertor Group
	double P_he_pump_recovered	 = P_pump_divertor;
						  
    double P_neutron_DivGroup_in = f_n_div			*P_neutron; // M applied below
    
    double P_condRad_DivGroup_in = f_alpha_div		*P_alpha_loss + // = 0
                                   f_rad_div		*P_rad_chamb  +
                                   f_rad_edge		*P_rad_edge   +
                                   P_cond_div;
    // same as "P_thermal_divertor" above except this has P_pump_divertor added
                                   
//cout << "P_neutron_DivGroup_in = " << P_neutron_DivGroup_in << "\n";
//cout << "P_condRad_DivGroup_in = " << P_condRad_DivGroup_in << "\n";


// Conduction and radiation powers are only deposited in the divertor plates.
     P_DivGroup[condRad][divPlates]  = P_condRad_DivGroup_in; // OK
     P_DivGroup[condRad][divBlanket] = 0.0;
     P_DivGroup[condRad][divShield]  = 0.0;
         
// 24-26, 54. Split the power from neutrons into fractions that go into divertor plates,
// blanket behind divertor and shield behind divertor. Include the effect of neutron 
// energy multiplication. 50%, 40%, 10%, (TBD 1%?)
    double f_div_plates;  fin_po.seekg(24*iskip_po + itxt_po +1); fin_po >> f_div_plates; 
    double f_div_blanket; fin_po.seekg(25*iskip_po + itxt_po +1); fin_po >> f_div_blanket; 
    double f_div_shield;  fin_po.seekg(26*iskip_po + itxt_po +1); fin_po >> f_div_shield;
	//double f_div_LT; 	  fin_po.seekg(54*iskip_po + itxt_po +1); fin_po >> f_div_LT;
    
    P_DivGroup[neutron][divPlates]  = f_div_plates  * f_nmul_divertor    * P_neutron_DivGroup_in;
    P_DivGroup[neutron][divBlanket] = f_div_blanket * f_nmul_div_blanket * P_neutron_DivGroup_in;
    P_DivGroup[neutron][divShield]  = f_div_shield  * f_nmul_div_shield  * P_neutron_DivGroup_in;
    // Add LT VV to this P_DivGroup??

// Low-temp components (e.g. vacuum vessel)
// Need M?
//	double LT_div_loss = f_div_LT*P_thermal_divertor*f_nmul_div_LT;


// ====================================================================================
// ====================================================================================
// Add up Thermal Power

      double P_thermal      = 0.0; // intialize here
      double P_thermal_fwbs = 0.0;
      double P_thermal_div  = 0.0;

// Add powers collected from the FWBS group (cooled by Pb-Li for SiC, dual-coolant for DCLL)
// SiC  has ~27% PbLi by volume in FW, ~83% PbLi by volume in blanket... 
// DCLL has ~65% He by volume in FW, ~77% PbLi and ~12% He by volume in blanket...
        for(int i=0; i<2; i++)
      { for(int j=0; j<2; j++)
      { for(int k=0; k<2; k++)
      { P_thermal      = P_thermal      + P_FWBSGroup[i][j][k]; 
        P_thermal_fwbs = P_thermal_fwbs + P_FWBSGroup[i][j][k]; 
//cout << i << j << k << P_FWBSGroup[i][j][k] << "\n";
} } }

   
// Add powers collected from the Divertor group (cooled by He only for both SCLL and DCLL):
// SCLL has ~53% He by volume
// DCLL has ~53% He by volume
        for(int i=0; i<2; i++)
      { for(int j=0; j<3; j++)
      { P_thermal     = P_thermal     + P_DivGroup[i][j]; // continue to add to P_thermal from above
        P_thermal_div = P_thermal_div + P_DivGroup[i][j]; } }

// C.K: Pthermal = 1.1*Pneutron + Palpha + Paux
//cout << "P_thermal CK = " << 1.1*P_neutron + P_alpha + 64.0;

// Add blanket and divertor pumping powers to P_thermal
// This is where the blanket pumping power is recovered as heat
	P_thermal_fwbs = P_thermal_fwbs + P_lm_pump_recovered;
	P_thermal_div  = P_thermal_div  + P_he_pump_recovered;
	P_thermal      = P_thermal_fwbs + P_thermal_div;

//cout << "P_thermal_fwbs   	= " << P_thermal_fwbs/1.0e6 << "\n";
//cout << "P_thermal_div    	= " << P_thermal_div/1.0e6 << "\n";
//cout << "P_thermal (total) 	= "	<< P_thermal/1.0e6 << "\n";
//cout << "P_thermal_div fraction	= " << P_thermal_div/P_thermal << "\n"; // ~28% DCLL, 25% SiC

// ====================================================================================
// ====================================================================================
// Gross Electric Power: DP ~920

	double eta_brayton      = brayton(flux_neutron_wall_Max, flux_heat_wall_Max, blanketmode);		 	 
	double P_gross_electric = eta_brayton*P_thermal;
	double P_brayton_loss   = (1.0 - eta_brayton)*P_thermal; 
	//which P_thermal? w/ P_pump_blanket thermal contribution and w/o "recovered" pump power


// ====================================================================================
// Thermal power recovered from pumps

//  1. LM MHD pump
//  ...not in use. Siegfried says that MHD pumps will never get 2 MPa head
/*	double P_lm_pump_recovered = P_pump_blanket/eta_pump_lm_mhd * 
								 eta_pump_lm_mhd_recovered      * eta_brayton; */

// NEED breakup of pumping ratios for DCLL dual-coolant...
//  2. LM mechanical pump
//	P_lm_pump_recovered = P_pump_blanket/eta_pump_lm_mech  * 
//						  eta_pump_lm_mech_recovered       * eta_brayton;
  
//  3. He pump/blower
//	P_he_pump_recovered = P_pump_divertor/eta_pump_he      * 
//						  eta_pump_he_recovered            * eta_brayton;

//cout << "P_lm_pump_recovered 	= " << P_lm_pump_recovered/1.0e6 << "\n";
//cout << "P_he_pump_recovered 	= " << P_he_pump_recovered/1.0e6 << "\n";

// Note - add Ohmic power someday.  But our plasma now don't have any.

// ====================================================================================
// ====================================================================================
// Net electric power:      (P_recirc = P_gross - P_net)
     double P_net_electric = P_gross_electric
      				   	   - P_plasma_heating/ eta_plasma_heating
                           - P_cd_generic/     eta_cd_generic
/*                         - P_cd[nb]/         eta_cd[nb]
                           - P_cd[lh]/         eta_cd[lh]
                           - P_cd[icrf]/       eta_cd[icrf]
                           - P_cd[ec]/         eta_cd[ec]
Commented out b/c the physics code doesn't have them yet. To be implemented by Chuck Kessel. */
                           - P_aux_func
                           - P_cryo
                           - P_pump_blanket/   eta_pump_lm_mech
                           - P_pump_divertor/  eta_pump_he;


// LCC might be good to print out wall-plug electrical power required with efficiencies
// Electrical power required to pump blanket = w/ eta, more than pp
// Blanket pumping power = w/o eta, given by eqn

//cout << "P_pump_blanket 		= " << P_pump_blanket << "\n";
//cout << "P_lm_pump_recovered 	= " << P_lm_pump_recovered << "\n";
//cout << "P_pump_divertor 		= " << P_pump_divertor << "\n";
//cout << "P_he_pump_recovered 	= " << P_he_pump_recovered << "\n";
                         

// ====================================================================================
// Power Losses Summary:

// P_fusion*(0.8*M + 0.2) = P_net_electric + P_losses (include P_aux and P_cryo?)
// Should be 8 losses, Brayton is largest
	double P_generation = P_fusion*(0.8*1.1 + 0.2)/1.0e6;
	double P_losses = P_plasma_heating_loss
					+ P_cd_loss
					+ P_brayton_loss
					+ P_pump_blanket_loss
					+ P_pump_divertor_loss
					+ P_aux_func
					+ P_cryo;
//					+ LT_fwbs_loss
//					+ LT_div_loss; // in W

	double P_leftover = P_generation - (P_net_electric/1.0e6) - (P_losses/1.0e6);

// Raise a flag if the power flows do not balance
	if(P_leftover > 2.0  || P_leftover < -2.0)
	{ cout << "Power flows do not balance! P_leftover = " << P_leftover << "\n"; }

/*cout << " P_plasma_heating_loss 	= " << P_plasma_heating_loss/1.0e6 << "\n";
cout << " P_cd_loss		= " << P_cd_loss/1.0e6 << "\n";
cout << " P_brayton_loss 	= " << P_brayton_loss/1.0e6 << "\n";
cout << " P_pump_blanket_loss 	= " << P_pump_blanket_loss/1.0e6 << "\n";
cout << " P_pump_divertor_loss 	= " << P_pump_divertor_loss/1.0e6 << "\n";
cout << " P_aux_func+P_cryo	= " << (P_aux_func+P_cryo)/1.0e6 << "\n";
//cout << " LT_fwbs_loss 		= " << LT_fwbs_loss/1.0e6 << "\n";
//cout << " LT_div_loss 		= " << LT_div_loss/1.0e6 << "\n";
cout << "  P_losses 		= " << P_losses/1.0e6 << "\n";
cout << "  P_net_elec 		= " << P_net_electric/1.0e6 << "\n";
cout << "  P_generation		 = " << P_generation << "\n";
cout << "   Leftover 		= " <<  P_generation-(P_losses/1.0e6)-(P_net_electric/1.0e6) << "\n";
*/


// P_nelec filter for 1 GW +/- 100 MW. Just puts a zero at the 7th filter location in putput file.
      if(P_net_electric > 1.1*1.e9 || P_net_electric <  0.9*1.e9) get_data_point_filter(7, 7);
      
      
// LCC does this below need to be turned on??  Or are we solely using Max heat flux, not Ave?
// Average Heat Flux on Divertor, does it need *2 to mirror top/bottom?
/* 	double A_divertor = 16.0*(*Div)[0].surface_area("outer") +
          	        	16.0*(*Div)[1].surface_area("outer") +
                      	16.0*(*Div)[2].surface_area("outer");
    double flux_divertor_Ave = P_particle/A_divertor; 
*/

// ====================================================================================
// FPC Volume Calulated
// used in DP.cpp ~1560, CA.cpp ~807

// Calc divertor volume
  double DivVol = 0.0; 
  	for(int i=0; i<6; i++)
  		{ DivVol = DivVol + (*Div)[i].volume(); }
  	
// Calc PF coils + spares including central solenoid
  double PoloFCVolume = 0.0; double PoloFCVolumeSpares = 0.0;
  	for (int i=0; i<36; i++)
  		{ PoloFCVolume = PoloFCVolume + (*PoloFCoil)[i].volume(); }  
  	for(int i=4; i<18; i++) 
  		{ PoloFCVolumeSpares = PoloFCVolumeSpares + (*PoloFCoil)[i].volume(); }

// Penetration shielding thickness (also read in CA ~338, DP ~1466)
// NEW - Read in penetration shielding thickness from powerflow.data input file
  double pen_shielding; fin_po.seekg(43*iskip_po + itxt_po +1); fin_po >> pen_shielding;
  double fpc_vol;
  fin_po.close();
  
if (blanketmode == "dcll") { 
	 	 fpc_vol = ((*FirstWallInb).volume() +
                    (*FirstWallOutb).volume() +
                    (*BackWallInboard).volume() +
                    (*BackWallOutboard).volume() +
                    (*BlanketInboard).volume() +
                    (*BlanketOutboard).volume() +
                    (*RepHTShield)[0].volume() +
                    (*RepHTShield)[1].volume() +
                    (1.0+pen_shielding)*((*HTShieldDCLL).volume()) +
                    (*SkeletonRing).volume() +                    
                    (*VacuumVesselDCLL)[0].volume() + 
                    (*VacuumVesselDCLL)[1].volume() +
                    (*VacuumVesselDCLL)[2].volume() + 
                    (*VacuumVesselDCLL)[3].volume() +
                    DivVol + (*TFCoil).volume() +
                    PoloFCVolume + PoloFCVolumeSpares)*16.0; } // everything *16

else {	 fpc_vol = ((*FirstWallInb).volume() +
                    (*FirstWallOutb).volume() +
                    (*BlanketInboard).volume() +
                    (*BlanketOutboard).volume() +
                    (*BlanketII).volume() +
                    (*RepHTShield)[0].volume() +
                    (*RepHTShield)[1].volume() +
                    (1.0+pen_shielding)*((*HTShield)[0].volume() + (*HTShield)[1].volume()) +
                    (*VacuumVessel)[0].volume() + 
                    (*VacuumVessel)[1].volume() +
                    DivVol + (*TFCoil).volume() +
                    PoloFCVolume + PoloFCVolumeSpares)*16.0; } // everything *16
                    
	get_fpc_vol(fpc_vol); // puts fpc_vol into this fcn for calling elsewhere
//cout << "FPC volume1 = " << fpc_vol << "\n"; // ~1220 m3


// *************************************************************************************
// LCC max heat flux on divertor moved from here to line 338

// ====================================================================================
// Inputs to functions for use in CA.cpp
	get_gross_electrical_power(P_gross_electric/1.e6);
	get_net_electrical_power  (P_net_electric/1.e6);
	get_thermal_power         (P_thermal/1.e6);
	get_fusion_power          (P_fusion/1.e6);
	get_P_pump_blanket        (P_pump_blanket/1.e6);
	get_P_pump_blanket        (P_pump_blanket/1.e6);
	
	
// NEW - Read in divertor, blanket design from setup.data input file
	ifstream fin_div;  fin_div.open("InputData/setup.data");
	fin_div.width(10); fin_div.precision(2);
	fin_div.setf(ios_base::uppercase, ios_base::scientific);
	int itxt_div = 28;  // length of line
	int toggle_ttube;  fin_div.seekg(5*itxt_div); fin_div >> toggle_ttube;
	int toggle_plate;  fin_div.seekg(7*itxt_div); fin_div >> toggle_plate;
	int toggle_finger; fin_div.seekg(9*itxt_div); fin_div >> toggle_finger;

	fin_div.close();

// ====================================================================================
// Printout when "system" mode is selected Powers.data
if(dispmode == "system") 
   { string filedir, filepath, extension; filedir = "./OutputParts/"; extension = ".data";
   filepath = filedir + "Powers" + extension; ofstream fileout(filepath.c_str()); 
   time_t rawtime;	struct tm * timeinfo; time ( &rawtime );
   timeinfo = localtime ( &rawtime );

// Output format is revised to allow easy alignment in Excel using "tab;" as delimiters
 fileout << "Powers [MW]: \n";
 fileout << "Printed on " 							<< asctime (timeinfo) << "\n";
 fileout << " Plasma Power [MW]:	;" 				<< P_plasma/1.0e6 << "\n";
 fileout << " Alpha Power [MW]:	;" 					<< P_alpha/1.0e6 << "\n";
 fileout << " Neutron Power [MW]:	;" 				<< P_neutron/1.0e6 << "\n";
 fileout << " Thermal Power [MW]:	;" 				<< P_thermal/1.0e6 << "\n";
 fileout << "  Thermal Power - FWBS:	;" 			<< P_thermal_fwbs/1.0e6 << "\n";
 fileout << "  Thermal Power - divertor:	;" 		<< P_thermal_div/1.0e6 << "\n";
 fileout << " Fusion Power [MW]:	;" 				<< P_fusion/1.0e6 << "\n";
 fileout << " Gross Electrical Power [MW]:	;" 		<< P_gross_electric/1.0e6 << "\n";
 fileout << " Net Electrical Power [MW]:	;" 		<< P_net_electric/1.0e6 << "\n";
 fileout << "\n";

// Power Balance
 fileout << "Power Balance Check: \n";
 fileout << " Power Generation [MW]:	;"			<< P_generation << "\n";
 fileout << " P_losses [MW]: \n";
 fileout << "  P_plasma_heating_loss [MW]:	;"		<< P_plasma_heating_loss/1.0e6 << "\n";
 fileout << "  P_cd_loss [MW]:	;" 					<< P_cd_loss/1.0e6 << "\n";
 fileout << "  P_brayton_loss [MW]:	;" 				<< P_brayton_loss/1.0e6 << "\n";
 fileout << "  P_pump_blanket_loss [MW]:	;"		<< P_pump_blanket_loss/1.0e6 << "\n";
 fileout << "  P_pump_divertor_loss [MW]	;"		<< P_pump_divertor_loss/1.0e6 << "\n";
 fileout << "  P_aux_func+P_cryo [MW]:	;"			<< (P_aux_func+P_cryo)/1.0e6 << "\n";
 //fileout << "  LT_fwbs_loss [MW]:	;"				<< LT_fwbs_loss/1.0e6 << "\n";
 //fileout << "  LT_div_loss [MW]:	;"				<< LT_div_loss/1.0e6 << "\n";
 fileout << " Sum of P_losses [MW]:	;" 				<< P_losses/1.0e6 << "\n";
 fileout << " Net Electrical Power [MW]:	;" 		<< P_net_electric/1.0e6 << "\n";
 fileout << " Leftover (should = 0) [MW]:	;" 		<< P_leftover << "\n";
 fileout << "\n";




// Return divertor in use
 fileout << "Divertor: \n";
if      (toggle_ttube  == 1) { fileout << " Divertor in use:	;" << "T-tube design, Dec. 2010, helium-cooled" << "\n"; }
else if (toggle_plate  == 1) { fileout << " Divertor in use:	;" << "plate-type (HCFC), helium-cooled" << "\n"; }
else if (toggle_finger == 1) { fileout << " Divertor in use:	;" << "combined plate-finger-type (HCPF), helium-cooled" << "\n"; }
 fileout << " He pump/blower efficiency:	;" 			<< eta_pump_he << "\n";
 fileout << " Divertor pumping power (He only):	;"		<< P_pump_divertor/1.0e6 << "\n";
 fileout << " Divertor pumping power (electrical) (He only):	;"	<< P_pump_divertor/eta_pump_he/1.0e6 << "\n";
 fileout << " Recovered power from He pump [MW]:	;"	<< P_he_pump_recovered/1.0e6 << "\n";
 //fileout << "  (P_he_pump_recovered = P_pump_divertor / eta_pump_he * eta_pump_he_recovered * eta_brayton)\n\n";
 fileout << " Power to Divertor [MW]:	;"				<< P_div << "\n";  
 fileout << " Divertor Heat Flux max [MW/m2]:	;"		<< flux_divertor_Max/1.0e6 << "\n"; 
 fileout << "\n";


// Return blanket in use
 fileout << "Blanket: \n";
if(blanketmode == "sic")  { fileout << " Blanket in use:	;" << "SCLL SiC blanket, PbLi-cooled, Aries-AT design" << "\n";
 fileout << " Blanket Version:	;" 					<< "** to do: add blanket version designator **" << "\n";
 fileout << " LM pump in use:	;" 					<< "mechanical LM pump" << "\n"; // line ~160, 677
 fileout << " LM pump efficiency:	;" 				<< eta_pump_lm_mech << "\n";
 fileout << " Blanket pumping power (LM only):	;"	<< P_pump_blanket/1.0e6 << "\n";
 fileout << " Blanket pumping power (electrical) (LM only):	;"	<< P_pump_blanket/eta_pump_lm_mech/1.0e6 << "\n";
 fileout << " Recovered power from LM pump:	;"  	<< P_lm_pump_recovered/1.0e6 << "\n";
 //fileout << "  (P_lm_pump_recovered = P_pump_blanket / eta_pump_lm_mech * eta_pump_lm_mech_recovered * eta_brayton)\n";
 fileout << "\n";
}

if(blanketmode == "dcll") { fileout << " Blanket in use:	;" << "DCLL blanket, PbLi and He-cooled, April 2009 design" << "\n";
 fileout << " Blanket Version:	;" 					<< "** to do: add blanket version designator **" << "\n";
 fileout << " LM pump in use:	;" 					<< "mechanical LM pump" << "\n"; // line ~160, 677
 fileout << " LM pump efficiency:	;" 				<< eta_pump_lm_mech << "\n";
 fileout << " He pump/blower efficiency:	;" 		<< eta_pump_he << "\n";
 fileout << " Blanket pumping power (LM+He combined):	;"	<< P_pump_blanket/1.0e6 << "\n";
 fileout << " Blanket pumping power (electrical) (LM+He combined):	;"	<< P_pump_blanket/eta_pump_lm_mech/1.0e6 << "\n";
 fileout << " Recovered power from both LM+He pumps:	;"  << P_lm_pump_recovered/1.0e6 << "\n";
 // Need breakup of PbLi/He pumping ratios, can't include both yet
 //fileout << " Recovered power from He pump:	;" 		<< P_he_pump_recovered/1.0e6 << "\n";
 //fileout << "  (P_lm_pump_recovered = P_pump_blanket / eta_pump_lm_mech * eta_pump_lm_mech_recovered * eta_brayton)\n"; 
 fileout << "\n";
}

 fileout << "Magnets: \n";
 fileout << " TF magnets in use:	;" 				<< "low-temp SC magnet technology" << "\n"; 
 fileout << " Nb3Sn at 4.2 K, 'Next Step, conservative costing' \n\n"; // from Geometry ~1850
 
 fileout << "Additional Parameters: \n";
 fileout << " Gross Cycle Efficiency (Brayton):	;"	<< eta_brayton << "\n";
 fileout << " Plasma Heating [MW]:	;"	 			<< P_plasma_heating/1.0e6 << "\n";
 fileout << " Plasma Heating Efficiency:	;" 		<< eta_plasma_heating << "\n";
 fileout << " Plasma Heating (electrical) [MW]:	;"	<< P_plasma_heating/eta_plasma_heating/1.0e6 << "\n";
 fileout << " Current Drive Power [MW]:	;" 			<< P_cd_generic/1.0e6 << "\n";
 fileout << " Current Drive Efficiency:	;" 			<< eta_cd_generic << "\n";
 fileout << " Current Drive Power (electrical) [MW]:	;" << P_cd_generic/eta_cd_generic/1.0e6 << "\n";
 fileout << " Aux. Functions Power [MW]:	;"		<< P_aux_func/1.0e6 << "\n";
 fileout << " Cryo Power [MW]:	;"	 				<< P_cryo/1.0e6 << "\n\n";
 
 fileout << " NWL ave (plasma surface) [MW/m2]:	;"	<< flux_neutron_wall_Ave/1.0e6 << "\n";
 fileout << " NWL max (plasma surface) [MW/m2]:	;"	<< flux_neutron_wall_Max/1.0e6 << "\n";
 fileout << " NWL ave (incl. SOL) [MW/m2]:	;"		<< flux_neutron_wall_Ave_w_SOL/1.0e6 << "\n";
 fileout << " NWL max (incl. SOL) [MW/m2]:	;"		<< flux_neutron_wall_Ave_w_SOL*peaking_factor/1.0e6 << "\n";
 fileout << " Peaking factor:	;"					<< peaking_factor << "\n\n";
  
 fileout << " Power to FW [MW]:	;"					<< P_rad_first_wall/1.0e6 << "\n";
 fileout << " FW Heat Flux ave [MW/m2]:	;"			<< flux_heat_wall_Ave/1.0e6 << "\n";
 fileout << " FW Heat Flux max [MW/m2]:	;"			<< flux_heat_wall_Max/1.0e6 << "\n\n";
} }

// END of "get_power_flow" function


// ====================================================================================
// ====================================================================================
// GROSS CYCLE EFFICIENCY of Brayton cycle
// for SiC Blanket "single coolant SiC/PbLi" and DCLL Blanket "dual-coolant PbLi"

// Reference: R. Raffray of MAX neutron wall load (NWL) at first wall vs 
// gross cycle efficiency (Q) with different max surface heat flux q"max; 
// See "Raffray - Pumping Power systems code input.xls"
// LCC removed old tables 5/9/11 and placed in "removedcode.cpp"

// 2010-11 M. Tillack concluded that we can continue to use R. Raffray's NWL vs cycle efficiency
// for the near term since our NWLmax and q" are low (~4.2 MW/m^2, ~0.25 MW/m^2 for ~AT).

double DesignPoint::brayton(double nwl, double q, string blanketmode)
{	double eta;
	double q2   = q/1.0e6; // change q, nwl to MW/m2
	double nwl2 = nwl/1.0e6;

if(blanketmode == "sic") {

// Polynomial fit to Raffray's Brayton table for SiC blanket
if (q2 <  0.250) { eta = 0; cout << "Maximum surface heat flux q'' is too low for the Brayton cycle tabulated interval " << q2 << "\n"; }

if (q2 >= 0.250 && q2 <= 0.375) { eta = -0.00347*nwl2 + 0.58859; 
if (nwl2 < 1.5 && nwl2 > 6.0) cout << "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " << nwl2 << "\n"; }

if (q2 >  0.375 && q2 <= 0.500) { eta = -0.00361*nwl2 + 0.56934; 
if (nwl2 < 1.8 && nwl2 > 6.0) cout << "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " << nwl2 << "\n"; }

if (q2 >  0.500 && q2 <= 0.625) { eta = -0.00362*nwl2 + 0.54828; 
if (nwl2 < 2.4 && nwl2 > 6.0) cout << "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " << nwl2 << "\n"; }

if (q2 >  0.625 && q2 <= 0.750) { eta = -0.00363*nwl2 + 0.52605; 
if (nwl2 < 3.0 && nwl2 > 6.0) cout << "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " << nwl2 << "\n"; }

if (q2 >  0.750) { eta = 0; cout << "Maximum surface heat flux q'' (" << q2 << ") is too high for the Brayton cycle tabulated interval.\n";}

//cout << "q''max = " << q2 << "\n";
//cout << "NWL_max = " << nwl2 << "\n";
//cout << "Brayton eta sic = " << eta << "\n";

return eta; 

} // end of SiC


// ***********OR**************OR*****************OR***************OR*******************

// Gross Cycle Efficiency for DCLL Blanket "dual-coolant PbLi" - Brayton cycle

if(blanketmode == "dcll") {

// Polynomial fit to Raffray's Brayton table for DCLL blanket
if (q2 <  0.250) { eta = 0; cout << "Maximum surface heat flux q'' is too low for the Brayton cycle tabulated interval " << q2 << "\n"; }

if (q2 >= 0.250 && q2 <= 0.500) { eta = 0.00013*pow(nwl2,3) - 0.00190*pow(nwl2,2) + 0.00087*nwl2 + 0.44578; 
if (nwl2 < 2.0 && nwl2 > 7.5) cout << "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " << nwl2 << "\n"; }

if (q2 >  0.500 && q2 <= 0.750) { eta = 0.00002*pow(nwl2,3) - 0.00009*pow(nwl2,2) - 0.00814*nwl2 + 0.45770; 
if (nwl2 < 2.4 && nwl2 > 7.5) cout << "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " << nwl2 << "\n"; }

if (q2 >  0.750 && q2 <= 1.000) { eta = 0.00018*pow(nwl2,3) - 0.00275*pow(nwl2,2) + 0.00599*nwl2 + 0.43195; 
if (nwl2 < 4.5 && nwl2 > 7.5) cout << "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " << nwl2 << "\n"; }

if (q2 >  1.000 && q2 <= 1.250) { eta = -0.00031*pow(nwl2,2) + 0.00290*nwl2 + 0.38929; 
if (nwl2 < 4.8 && nwl2 > 7.5) cout << "Maximum neutron wall load (NWL_max) at first wall is out of the Brayton cycle tabulated interval " << nwl2 << "\n"; }

if (q2 >  0.750) { eta = 0; cout << "Maximum surface heat flux q'' (" << q2 << ") is too high for the Brayton cycle tabulated interval.\n";}

//cout << "q''max = " << q2 << "\n";
//cout << "NWL_max = " << nwl2 << "\n";
//cout << "Brayton eta dcll = " << eta << "\n";

return eta; 

} // end of DCLL
} // end of Brayton fcn
	


// ****** NEW function *****========================================================
// =================================================================================
// Divertor Pumping Power Correlations

// All divertors are He-cooled; function is called by line ~675
// Reference: "Divertor Pumping Powers - ppcompare2.xls" Dec. 2010 by M. Tillack, X. Wang

double DesignPoint::pumping_power_divertor(double flux_divertor_Max)
{
// 1.) J. BURKE'S T-Tube 2010 divertor design pumping power curve
// Linear-taper with fixed slot width of 0.45 mm, inlet temp 600C, outlet 677C, P = 10 MPa
// surface heat flux (MW/m^2) vs P_pump/P_thermal (%)
// NEED div max heat flux, "flux_divertor_Max" = essentially "Qdivoutb"
	double ppttube = 0.05594*pow(flux_divertor_Max/1.0e6,2) + 
					 0.04931*(flux_divertor_Max/1.0e6) - 0.08702;
	if (ppttube<0.0)
	    ppttube = 0.0;

// 2.) X. WANG'S plate divertor design (HCFP) *most common*
// Helium-cooled flat plate, inlet temp 600C, outlet 677C, P = 10 MPa
	double ppplate = 0.0392*pow(flux_divertor_Max/1.0e6,3) - 
					 0.743 *pow(flux_divertor_Max/1.0e6,2) +
					 5.44  *(flux_divertor_Max/1.0e6) - 11.1;
	if (ppplate<0.0)
	    ppplate = 0.0;
	 			
// 3.) X. WANG'S finger divertor design	(HCPF)
// Helium-cooled combined plate and finger, inlet temp 600C, outlet 700C, P = 10 MPa
	double ppfinger = 0.0335*pow(flux_divertor_Max/1.0e6,2) + 
				 	  0.22  *(flux_divertor_Max/1.0e6) - 0.892;
	if (ppfinger<0.0)
	    ppfinger = 0.0;

// 4.) M. YODA'S Georgia Tech Divertor Correlation to HCPF w/ and w/out fins
// HCFP correlations from Yoda's July 2011 presentation
	double ppfingerGTbare = 8.91e-5*exp(0.502*flux_divertor_Max/1.0e6); // bare, 2 mm slots
//	double ppfingerGTfins = 1.12e-10*exp(1.102*flux_divertor_Max/1.0e6);// fins, 2 mm slots (not realistic)
	
//cout << "flux_divertor_Max/1.0e6  = " << flux_divertor_Max/1.0e6 << "\n";
//cout << "Pumping power T-tube div = " << ppttube << "\n";
//cout << "Pumping power plate div  = " << ppplate << "\n";
//cout << "Pumping power finger div = " << ppfinger << "\n";
//cout << "Pumping power GT bare    = " << ppfingerGTbare << "\n";

// Switch for different divertor designs, 1/0 to turn on/off
// NEW - Read in divertor design from setup.data input file
	ifstream fin_div;  fin_div.open("InputData/setup.data");
	fin_div.width(10); fin_div.precision(2);
	fin_div.setf(ios_base::uppercase, ios_base::scientific);
	int itxt_div = 28;  // length of line
	
	int toggle_ttube;  fin_div.seekg(5*itxt_div); fin_div >> toggle_ttube;
	int toggle_plate;  fin_div.seekg(7*itxt_div); fin_div >> toggle_plate;
	int toggle_finger; fin_div.seekg(9*itxt_div); fin_div >> toggle_finger;

	fin_div.close();

// Return pumping power, output is the % of P_thermal of He ONLY
	if      (toggle_ttube  == 1) { return ppttube; }
	else if (toggle_plate  == 1) { return ppplate; }
	else if (toggle_finger == 1) { return ppfinger; }

// Correlations are valid from q'' 5-15 MW/m2, filter applied above when eqn is called

} // eof


// Old SiC pumping power tables were here, removed to "removedcode" 5/9/11 LCC
// ====================================================================================
// ====================================================================================

// FPC_volume transfer to CA.cpp, LCC added
void DesignPoint::get_fpc_vol(double volume_fpc)
	{ fpcvolume = volume_fpc; }

void DesignPoint::get_gross_electrical_power(double power)
	{ power_gross_electrical = power; }

void DesignPoint::get_net_electrical_power(double power)
	{ power_net_electrical = power; }

void DesignPoint::get_thermal_power(double power)
	{ power_thermal = power; }

void DesignPoint::get_P_pump_blanket(double power)
	{ power_blanket = power; }

void DesignPoint::get_fusion_power(double power)
	{ power_fusion = power; }

void DesignPoint::get_physics(double nwl, double qmhd, double aip, double betap)
	{ ave_neutron_wall_load = nwl;
	  q95 					= qmhd;
 	  plsm_current 			= aip;
 	  plsm_poloidal_beta 	= betap; }
	
void DesignPoint::get_neutron_wall_loads(double average, double maximum)
	{ ave_neutron_wall_load = average;
   	  max_neutron_wall_load = maximum; }

void DesignPoint::get_divertor_loads(double inboard, double outboard)
	{ max_heat_flux_div_outboard = outboard;
      max_heat_flux_div_inboard  = inboard; }

void DesignPoint::get_fw_heat_fluxes(double average, double maximum)
	{ ave_heat_flux_fw = average;
   	  max_heat_flux_fw = maximum; }

void DesignPoint::get_data_point_filter(int dpf, int loc)
	{ data_flag[loc] = dpf;   // dpf is the data point filter, should be applied to loc > 0
   	  data_flag[0]   = dpf; } // if loc is 0, assign the latest filter value to it


// ============================================================================================
// ============================================================================================
// Detailed System.data printout when "SYSTEM" mode is selected, called by Aries.cpp

void DesignPoint::show
			(Part * Plasma, 			Part * Separatrix,			Part * FWInboard,
			 Part * FWOutboard,			Part * BlanketInboard,		Part * BlanketOutboard,
			 Part * BlanketII, 			Part (*Divertor)[6],		Part (*RepHTShield)[2],
			 Part (*HTShield)[2],		Part (*VVessel)[4],			Part * TFC,
			 Part * Window,				Part * BCyl,				Part * UpperCap,
			 Part * LowerCap,			Part (*PoloFC)[36],			Part * Cryo, 
             Part * MPort, 				Part (*Duct)[3],			Part * Plasma_AT_w_SOL,
             double (*Outputs)[100],	Part * Plasma_DCLL,			Part * SkeletonRing,
             Part * BackWallInboard,	Part * BackWallOutboard,	Part (*VacuumVesselDCLL)[4],
             Part * HTShieldDCLL,		Part * Manifold,			Part (*Manifolding)[2])	
             // added for DCLL: Plasma_DCLL - HTShieldDCLL
{
//cout << plasma_params[3] << "\n"; // not working because show was using designpoint AT not DCLL

// Display the Contours
// calls Part.cpp line ~1000 to create geometry files for Matlab visualization

// DCLL-only parts:
if(blanketmode == "dcll") {	
	(*Plasma_DCLL).show(); // bus error when running sic mode
	(*BackWallInboard).show();
	(*BackWallOutboard).show();
	for(int i=0; i<4; i++)(*VacuumVesselDCLL)[i].show();
	(*SkeletonRing).show();
	(*HTShieldDCLL).show();
	(*Manifold).show();
}

// SCLL-only parts:
else {	
	(*Plasma).show();
	//(*Plasma_AT_w_SOL).show(); // doesn't produce part	
	(*BlanketII).show();
	(*HTShield)[0].show();
	(*HTShield)[1].show();	
	(*VVessel)[0].show();
	(*VVessel)[1].show(); // show [2],[3]??
	(*Window).show();
	(*Manifolding)[0].show(); // split up here, each their own show()
	(*Manifolding)[1].show();
}
	
// Shared parts:
	(*Separatrix).show();
	(*FWInboard).show();
	(*FWOutboard).show();
	(*BlanketInboard).show();
	(*BlanketOutboard).show();
	for(int i=0; i<6; i++) (*Divertor)[i].show();
	(*RepHTShield)[0].show();
	(*RepHTShield)[1].show();
	(*TFC).show();
	(*BCyl).show();
	(*UpperCap).show();
	(*LowerCap).show();
	for(int i=0; i<36; i++) (*PoloFC)[i].show();
	(*Cryo).show();
	(*MPort).show();
	(*Duct)[0].show();
	(*Duct)[1].show();
	//(*Duct)[2].show(); // problematic, causes bus error for sic, perhaps there are no contours to show
	// plus the duct2 IBL output file is empty
	// BUS ERROR occurs when it tries to read something that isn't there in memory...	
    const double pi = 4.0*atan(1.0);


// Setup output file for mart compositions
   string filedir, filepath, extension;
   filedir = "./OutputParts/"; extension = ".data";
   filepath = filedir + "System" + extension;
   ofstream fileout(filepath.c_str());
   
// Printout date/time
	time_t rawtime;	struct tm * timeinfo; time ( &rawtime );
	timeinfo = localtime ( &rawtime );

// try fileout << setw(12);  setprecision(5), cout.width(10), fileout.width(15)


// =====================================================================
 fileout << "Summary of Power Plant System Parameters: \n";
 if(blanketmode == "dcll")  fileout << "-- DCLL Blanket -- \n";
 else 						fileout << "-- SCLL Blanket -- \n";

 fileout << "SVN Rev # (add manually),  Printed on " 					<< asctime (timeinfo) << "\n"; // make rev # an input somehow
 fileout << " Plasma major radius [m]:	;" 								<< (*Outputs)[0] << "\n";
 fileout << " Plasma aspect ratio [-]:	;" 								<< (*Outputs)[1] << "\n";
 fileout << " Toroidal magnetic field at plasma major radius [T]:	;" 	<< (*Outputs)[2] << "\n";
 fileout << " Plasma current [MA]:	;" 									<< (*Outputs)[3] << "\n";
 fileout << " Normalized Beta [-]:	;" 									<< (*Outputs)[4]*100.0 << "\n";
 fileout << " MHD safety factor [-]:	;" 								<< (*Outputs)[5] << "\n";
 fileout << " Cylindrical safety factor [-]:	;" 						<< (*Outputs)[6] << "\n";
 fileout << " Plasma triangularity [-]:	;" 								<< (*Outputs)[7] << "\n";
 fileout << " Exponent on plasma density profile [-]:	;" 				<< (*Outputs)[8] << "\n";
 fileout << " Exponent on plasma temperature profile [-]:	;" 			<< (*Outputs)[9] << "\n";
 fileout << " Greenwald density fraction [-]:	;" 						<< (*Outputs)[10] << "\n";
 fileout << " Fusion gain [-]:	;" 										<< (*Outputs)[11] << "\n";
 fileout << " Plasma elongation [-]:	;" 								<< (*Outputs)[12] << "\n";
 fileout << " H98(y,2) multiplier on energy confinement scaling [-]:	;" << (*Outputs)[13] << "\n";
 fileout << " Current diffusion time [s]:	;" 							<< (*Outputs)[14] << "\n";
 fileout << " Energy confinement time [s]:	;" 							<< (*Outputs)[15] << "\n";
 fileout << " Ratio of effec. particle conf. time to energy conf. time [-]:	;" << (*Outputs)[16] << "\n";
 fileout << " Volume average electron density [m^-3]:	;" 				<< (*Outputs)[17] << "\n";
 fileout << " Volume average temperature [KeV]:	;" 						<< (*Outputs)[18] << "\n";
 fileout << " Flattop time [s]:	;" 										<< (*Outputs)[19] << "\n";
 fileout << " L to H-mode threshold power [MW]:	;" 						<< (*Outputs)[20] << "\n";
 fileout << " Ave neutron wall load (at plasma surface) [MW/m^2]:	;"	<< (*Outputs)[21]/1.e6 << "\n";
 //fileout << " Ave neutron wall load (NWL_ave w/ SOL) [MW/m^2]:	;"		<< flux_neutron_wall_Ave_w_SOL << "\n";
 fileout << " Bremsstrahlung radiation loss [MW]:	;" 					<< (*Outputs)[22] << "\n";
 fileout << " Bootstrap current fraction [-]:	;" 						<< (*Outputs)[23] << "\n";
 fileout << " Current drive efficiency for generic heating/CD source [-]:	;" << (*Outputs)[24] << "\n";
 fileout << " Power from generic CD source (steady-state) [MW]:	;" 		<< (*Outputs)[25] << "\n";
 fileout << " Total aux. heating power into plasma for power balance [MW]:	;" << (*Outputs)[26] << "\n";
 fileout << " Cyclotron radiation loss [MW]:	;" 						<< (*Outputs)[27] << "\n";
 fileout << " Ohmic power in plasma [MW]:	;" 							<< (*Outputs)[28] << "\n";
 fileout << " Line radiation loss [MW]:	;" 								<< (*Outputs)[29] << "\n";
 fileout << " Fast alpha beta [-]:	;" 									<< (*Outputs)[30] << "\n";
 fileout << " Effective charge [-]:	;" 									<< (*Outputs)[31] << "\n";
 fileout << " Ratio of plasma temperature at center to edge [-]:	;" 	<< (*Outputs)[32] << "\n";
 fileout << " Ratio of plasma density at center to edge [-]:	;" 		<< (*Outputs)[33] << "\n";
 fileout << " Eta_CD for source 1 [-]:	;" 								<< (*Outputs)[34] << "\n";
 fileout << " Eta_CD for source 2 [-]:	;" 								<< (*Outputs)[35] << "\n";
 fileout << " Eta_CD for source 3 [-]:	;" 								<< (*Outputs)[36] << "\n";
 fileout << " Eta_CD for source 4 [-]:	;" 								<< (*Outputs)[37] << "\n";
 fileout << " Impurity fraction 1 [-]:	;" 								<< (*Outputs)[38] << "\n";
 fileout << " Impurity fraction 2 [-]:	;" 								<< (*Outputs)[39] << "\n";
 fileout << " Impurity fraction 3 [-]:	;" 								<< (*Outputs)[40] << "\n";
 fileout << " Peak temperature [keV]:	;" 								<< (*Outputs)[41] << "\n";
 fileout << " External CD fraction of total plasma current [-]:	;" 		<< (*Outputs)[42] << "\n";
 fileout << " Total non-inductive fraction of total plasma current [-]:	;"		<< (*Outputs)[43] << "\n";
 fileout << " Density fraction of helium (frac. of electron density) [-]:	;" 	<< (*Outputs)[44] << "\n";
 fileout << " Density fraction of DT fuel (frac. of electron density) [-]:	;" 	<< (*Outputs)[45] << "\n";
 fileout << " Plasma stored energy [MJ]:	;" 							<< (*Outputs)[46] << "\n";
 fileout << " Poloidal flux swing required to raise plasma current [?]:	;" 		<< (*Outputs)[47] << "\n";
 fileout << " Loop voltage in plasma [V]:	;" 							<< (*Outputs)[48] << "\n";
 fileout << " Density profile peak-to-volume average [-]:	;" 			<< (*Outputs)[49] << "\n";
 fileout << " Temperature profile peak-to-volume average [-]:	;" 		<< (*Outputs)[50] << "\n";
 fileout << " Poloidal beta [-]:	;" 									<< (*Outputs)[51] << "\n";
 fileout << " Toroidal beta [-]:	;" 									<< (*Outputs)[52] << "\n";
 fileout << " Alpha power in plasma [MW]:	;" 							<< (*Outputs)[53] << "\n";
 fileout << " Ratio of loss power to L-H threshold power [-]:	;" 		<< (*Outputs)[54] << "\n";
 fileout << " Plasma surface area [m^2]:	;"							<< (*Outputs)[55] << "\n"; // also A_plasma_surface
// fileout << " Plasma surface area (incl. SOL) [m^2]:	;"				<< A_plasma_surface_w_SOL << "\n";
 fileout << " Maximum magnetic field in the TF coil [T]:	;" 			<< (*Outputs)[56] << "\n";
 fileout << " Maximum IB heat flux on divertor [MW/m^2]:	;" 			<< (*Outputs)[57]/1.e6 << "\n";
 fileout << " Maximum OB heat flux on divertor [MW/m^2]:	;" 			<< (*Outputs)[58]/1.e6 << "\n";
 fileout << " Maximum heat flux to the first wall [MW/m^2]:	;" 			<< (*Outputs)[59]/1.e6 << "\n";
 fileout << " Gross electrical power [MW]:	;" 							<< (*Outputs)[60] + (*Outputs)[61] << "\n";
 fileout << " Net electric power [MW]:	;" 								<< (*Outputs)[60] << "\n";
 fileout << " Recirculated power (Pgross - Pnet) [MW]:	;" 				<< (*Outputs)[61] << "\n";
 fileout << " Thermal power [MW]:	;" 									<< (*Outputs)[62] << "\n";
 fileout << " Fusion power [MW]:	;" 									<< (*Outputs)[11]*(*Outputs)[26] << "\n";
 fileout << " Cost of electricity [mills/kWh]:	;" 						<< (*Outputs)[63] << "\n";
 fileout << " Data point number from plasma physics output file [-]:	;" << (*Outputs)[64] << "\n";
  fileout << "\n\n";


 fileout << "===============================================================\n";
 fileout << "Defining Parameters for Plasma Closed Surface Equation: \n\n";
 fileout << " Plasma major radius [m]:	;" 	<< plasma_params[0] << "\n";
 fileout << " Plasma minor radius [m]:	;" 	<< plasma_params[1] << "\n";
 fileout << " Elongation [-]:	;" 			<< plasma_params[2] << "\n";
 fileout << " Triangularity [-]:	;" 		<< plasma_params[3] << "\n";
 fileout << " Outer Squareness [-]:	;" 		<< plasma_params[4] << "\n";
 fileout << " Inner Squareness [-]:	;" 		<< plasma_params[5] << "\n";
 fileout << " Scrape-off layer [m]:	;" 		<< plasma_params[6] << "\n\n";


fileout << "Vertical Build [m]:\n"; 
// from Aries.cpp ~180, from plasma edge out to VV vertically
 fileout << " ... Work in progress, see 'builds' input file for the time being ...\n";
/* fileout << " SOL:	;" 						<< top_gen_bld[0] << "\n";
 fileout << " Divertor Plate:	;" 			<< top_gen_bld[0] << "\n";
 fileout << " Replaceable HT Shield:	;" 	<< top_gen_bld[1] << "\n";
 fileout << " Gap:	;" 						<< top_gen_bld[2] << "\n";
 fileout << " HT Shield:	;" 				<< top_gen_bld[3] << "\n";
if(blanketmode == "dcll") {
 fileout << " Skeleton Ring:	;" 			<< top_gen_bld[4] << "\n";
 fileout << " Gap:	;" 						<< top_gen_bld[5] << "\n";
 fileout << " Vacuum Vessel:	;" 			<< top_gen_bld[6] << "\n\n"; } 
else {
 fileout << " Gap:	;" 						<< top_gen_bld[4] << "\n";
 fileout << " Vacuum Vessel:	;" 			<< top_gen_bld[5] << "\n\n\n"; }
*/


fileout << "\n\n";
fileout << "===============================================================\n";
fileout << "Inboard to Outboard Radial Build [m]:\n";  
// NOTE: the "radius" is always the outermost dimension of the build

 fileout << "	;Thickness	;Radius \n";
 fileout << " Device Centerline	;"  << "0.0000" << "	;" << "0.0000" << "\n";
 fileout << " Bore	;" 				<< (*PoloFC)[0].bounding_box_equatorial(0, "outer") << "	;" << (*PoloFC)[0].bounding_box_equatorial(0, "outer") << "\n";
 fileout << " Central Solenoid	;"	<< (*PoloFC)[0].bounding_box_equatorial(1, "outer")-(*PoloFC)[0].bounding_box_equatorial(0, "outer") << "	;" 
                                    << (*PoloFC)[0].bounding_box_equatorial(1, "outer") << "\n";
 fileout << " Gap	;" 				<< (*BCyl).bounding_box_equatorial(0, "outer")-(*PoloFC)[0].bounding_box_equatorial(1, "outer") << "	;"
                                    << (*BCyl).bounding_box_equatorial(0, "outer") << "\n";
 fileout << " Bucking Cylinder	;" 	<< (*BCyl).bounding_box_equatorial(1, "outer")-(*BCyl).bounding_box_equatorial(0, "outer") << "	;"
                                    << (*BCyl).bounding_box_equatorial(1, "outer") << "\n";
 fileout << " Gap	;" 				<< (*TFC).bounding_box_equatorial(0, "outer")-(*BCyl).bounding_box_equatorial(1, "outer") << "	;"
                                    << (*TFC).bounding_box_equatorial(0, "outer") << "\n";
 fileout << " TF Coil Annulus	;" 	<< (*TFC).bounding_box_equatorial(0, "inner")-(*TFC).bounding_box_equatorial(0, "outer") << "	;"
                                    << (*TFC).bounding_box_equatorial(0, "inner") << "\n";
/*cout << (*PoloFC)[0].bounding_box_equatorial(0, "outer") << " PoloFC[0] 0, outer \n";
cout << (*PoloFC)[0].bounding_box_equatorial(1, "outer") << " PoloFC[0] 1, outer \n";
cout << (*BCyl).bounding_box_equatorial(0, "outer") << " BCyl 0, outer \n";
cout << (*BCyl).bounding_box_equatorial(1, "outer") << " BCyl 0, outer \n";
cout << (*TFC).bounding_box_equatorial(0, "outer") << " TFC 0, outer \n";
cout << (*TFC).bounding_box_equatorial(0, "inner") << " TFC 0, inner \n";
*/

if(blanketmode == "dcll") {
 fileout << " Gap	;" 				<< (*VacuumVesselDCLL)[0].bounding_box_equatorial(0, "outer")-(*TFC).bounding_box_equatorial(0, "inner") << "	;"
                                    << (*VacuumVesselDCLL)[0].bounding_box_equatorial(0, "outer") << "\n";
 fileout << " Vacuum Vessel IB	;"	<< (*VacuumVesselDCLL)[0].bounding_box_equatorial(1, "outer")-(*VacuumVesselDCLL)[0].bounding_box_equatorial(0, "outer") << "	;"
                                    << (*VacuumVesselDCLL)[0].bounding_box_equatorial(1, "outer") << "\n";
 fileout << " Gap	;"				<< (*SkeletonRing).bounding_box_equatorial(0, "outer")-(*VacuumVesselDCLL)[0].bounding_box_equatorial(1, "outer") << "	;"
									<< (*SkeletonRing).bounding_box_equatorial(0, "outer") << "\n";
 fileout << " Skeleton Ring IB	;"	<< (*SkeletonRing).bounding_box_equatorial(0, "inner")-(*SkeletonRing).bounding_box_equatorial(0, "outer") << "	;"
                                    << (*SkeletonRing).bounding_box_equatorial(0, "inner") << "\n";
 fileout << " Gap	;" 				<< (*HTShieldDCLL).bounding_box_equatorial(0, "outer")-(*SkeletonRing).bounding_box_equatorial(0, "inner") << "	;"
                                    << (*HTShieldDCLL).bounding_box_equatorial(0, "outer") << "\n";
 fileout << " HT Shield IB	;" 		<< (*HTShieldDCLL).bounding_box_equatorial(0, "inner")-(*HTShieldDCLL).bounding_box_equatorial(0, "outer") << "	;"
                                    << (*HTShieldDCLL).bounding_box_equatorial(0, "inner") << "\n"; 
 fileout << " Gap	;" 				<< (*BackWallInboard).bounding_box_equatorial(0, "outer")-(*HTShieldDCLL).bounding_box_equatorial(0, "inner") << "	;"
                                    << (*BackWallInboard).bounding_box_equatorial(0, "outer") << "\n";
 fileout << " Back Wall IB	;" 		<< (*BackWallInboard).bounding_box_equatorial(1, "outer")-(*BackWallInboard).bounding_box_equatorial(0, "outer") << "	;"
                                    << (*BackWallInboard).bounding_box_equatorial(1, "outer") << "\n";
 fileout << " Gap	;"				<< (*BlanketInboard).bounding_box_equatorial(0, "outer")-(*BackWallInboard).bounding_box_equatorial(1, "outer") << "	;"
                                    << (*BlanketInboard).bounding_box_equatorial(0, "outer") << "\n";
/*cout << (*VacuumVesselDCLL)[0].bounding_box_equatorial(0, "outer") << " VacuumVesselDCLL[0] 0, outer \n";
cout << (*VacuumVesselDCLL)[0].bounding_box_equatorial(1, "outer") << " VacuumVesselDCLL[0] 1, outer \n";
cout << (*SkeletonRing).bounding_box_equatorial(0, "outer") << " SkeletonRing 0, outer \n";
cout << (*SkeletonRing).bounding_box_equatorial(0, "inner") << " SkeletonRing 0, inner \n";
cout << (*HTShieldDCLL).bounding_box_equatorial(0, "outer") << " HTShieldDCLL 0, outer \n";
cout << (*HTShieldDCLL).bounding_box_equatorial(0, "inner") << " HTShieldDCLL 0, inner \n";
cout << (*BackWallInboard).bounding_box_equatorial(0, "outer") << " BackWallInboard 0, outer\n";
cout << (*BackWallInboard).bounding_box_equatorial(1, "outer") << " BackWallInboard 1, outer\n";
*/
}

else {
 fileout << " Gap	;" 				<< (*VVessel)[0].bounding_box_equatorial(0, "outer")-(*TFC).bounding_box_equatorial(0, "inner") << "	;"
                                    << (*VVessel)[0].bounding_box_equatorial(0, "outer") << "\n";
 fileout << " Vacuum Vessel IB	;"	<< (*VVessel)[0].bounding_box_equatorial(0, "inner")-(*VVessel)[0].bounding_box_equatorial(0, "outer") << "	;"
                                    << (*VVessel)[0].bounding_box_equatorial(0, "inner") << "\n"; 
 fileout << " Gap	;" 				<< (*HTShield)[0].bounding_box_equatorial(0, "outer")-(*VVessel)[0].bounding_box_equatorial(0, "inner") << "	;"
                                    << (*HTShield)[0].bounding_box_equatorial(0, "outer") << "\n";
 fileout << " HT Shield IB	;" 		<< (*HTShield)[0].bounding_box_equatorial(0, "inner")-(*HTShield)[0].bounding_box_equatorial(0, "outer") << "	;"
                                    << (*HTShield)[0].bounding_box_equatorial(0, "inner") << "\n";
 fileout << " Gap	;" 				<< (*BlanketInboard).bounding_box_equatorial(0, "outer")-(*HTShield)[0].bounding_box_equatorial(0, "inner") << "	;"
                                    << (*BlanketInboard).bounding_box_equatorial(0, "outer") << "\n";
/*cout << (*VVessel)[0].bounding_box_equatorial(0, "outer") << " VVessel[0] 0, outer \n";
cout << (*VVessel)[0].bounding_box_equatorial(0, "inner") << " VVessel[0] 0, inner \n";
cout << (*VVessel)[0].bounding_box_equatorial(1, "outer") << " VVessel[0] 1, outer \n";
cout << (*VVessel)[0].bounding_box_equatorial(1, "inner") << " VVessel[0] 1, inner \n";
cout << (*HTShield)[0].bounding_box_equatorial(0, "outer") << " HTShield[0] 0, outer \n";
cout << (*HTShield)[0].bounding_box_equatorial(0, "inner") << " HTShield[0] 0, inner \n";
cout << (*HTShield)[0].bounding_box_equatorial(1, "outer") << " HTShield[0] 1, outer \n";
cout << (*HTShield)[0].bounding_box_equatorial(1, "inner") << " HTShield[0] 1, inner \n";
cout << (*BlanketInboard).bounding_box_equatorial(0, "outer") << " BlanketInboard 0, outer \n";
cout << (*BlanketInboard).bounding_box_equatorial(1, "outer") << " BlanketInboard 1, outer \n";
cout << (*FWInboard).bounding_box_equatorial(0, "outer") << " FWInboard 0, outer \n";
cout << (*FWInboard).bounding_box_equatorial(1, "outer") << " FWInboard 1, outer \n";
*/
}

// Combined
 fileout << " Blanket IB	;" 		<< (*BlanketInboard).bounding_box_equatorial(1, "outer")-(*BlanketInboard).bounding_box_equatorial(0, "outer") << "	;"
                                    << (*BlanketInboard).bounding_box_equatorial(1, "outer") << "\n";
 fileout << " Gap	;" 				<< (*FWInboard).bounding_box_equatorial(0, "outer")-(*BlanketInboard).bounding_box_equatorial(1, "outer") << "	;"
                                    << (*FWInboard).bounding_box_equatorial(0, "outer") << "\n"; // LCC gap = 0
 fileout << " First Wall IB	;" 		<< (*FWInboard).bounding_box_equatorial(1, "outer")-(*FWInboard).bounding_box_equatorial(0, "outer") << "	;"
                                    << (*FWInboard).bounding_box_equatorial(1, "outer") << "\n";

double plasma_major_R;
if(blanketmode == "dcll") {
 fileout << " Scrape-Off Layer IB	;"	<< (*Plasma_DCLL).bounding_box(0, "outer")-(*FWInboard).bounding_box_equatorial(1, "outer") << "	;"
                                    	<< (*Plasma_DCLL).bounding_box(0, "outer") << "\n";
	plasma_major_R = ((*Plasma_DCLL).bounding_box(2, "outer") + (*Plasma_DCLL).bounding_box(0, "outer"))/2;
 fileout << " Plasma Centerline	;" 		<< ((*Plasma_DCLL).bounding_box(2, "outer") - (*Plasma_DCLL).bounding_box(0, "outer"))/2 << "	;"
                                    	<< plasma_major_R << "\n";
 fileout << " Plasma Outer Edge	;" 		<< ((*Plasma_DCLL).bounding_box(2, "outer") - (*Plasma_DCLL).bounding_box(0, "outer"))/2 << "	;"
                                    	<<  (*Plasma_DCLL).bounding_box(2, "outer") << "\n";
 fileout << " Scrape-Off Layer OB	;" 	<< (*FWOutboard).bounding_box_equatorial(0, "outer")-(*Plasma_DCLL).bounding_box(2, "outer") << "	;"
                                    	<< (*FWOutboard).bounding_box_equatorial(0, "outer") << "\n";
}

else {
 fileout << " Scrape-Off Layer IB	;"	<< (*Plasma).bounding_box(0, "outer")-(*FWInboard).bounding_box_equatorial(1, "outer") << "	;"
                                    	<< (*Plasma).bounding_box(0, "outer") << "\n";
	plasma_major_R = ((*Plasma).bounding_box(2, "outer") + (*Plasma).bounding_box(0, "outer"))/2;
 fileout << " Plasma Centerline	;" 		<< ((*Plasma).bounding_box(2, "outer") - (*Plasma).bounding_box(0, "outer"))/2 << "	;"
                                    	<< plasma_major_R << "\n";
 fileout << " Plasma Outer Edge	;" 		<< ((*Plasma).bounding_box(2, "outer") - (*Plasma).bounding_box(0, "outer"))/2 << "	;"
                                    	<<  (*Plasma).bounding_box(2, "outer") << "\n";
 fileout << " Scrape-Off Layer OB	;" 	<< (*FWOutboard).bounding_box_equatorial(0, "outer")-(*Plasma).bounding_box(2, "outer") << "	;"
                                    	<< (*FWOutboard).bounding_box_equatorial(0, "outer") << "\n";
 }
 
// Combined
 fileout << " First Wall OB	;" 		<< (*FWOutboard).bounding_box_equatorial(1, "outer")-(*FWOutboard).bounding_box_equatorial(0, "outer") << "	;"
                                    << (*FWOutboard).bounding_box_equatorial(1, "outer") << "\n";
 fileout << " Gap	;" 				<< (*BlanketOutboard).bounding_box_equatorial(0, "outer")-(*FWOutboard).bounding_box_equatorial(1, "outer") << "	;"
                                    << (*BlanketOutboard).bounding_box_equatorial(0, "outer") << "\n";
 fileout << " Blanket OB	;" 		<< (*BlanketOutboard).bounding_box_equatorial(1, "outer")-(*BlanketOutboard).bounding_box_equatorial(0, "outer") << "	;"
                                    << (*BlanketOutboard).bounding_box_equatorial(1, "outer") << "\n";
/*cout << (*FWOutboard).bounding_box_equatorial(0, "outer") << " FWOutboard 0, outer \n";
cout << (*FWOutboard).bounding_box_equatorial(1, "outer") << " FWOutboard 1, outer \n";
cout << (*BlanketOutboard).bounding_box_equatorial(0, "outer") << " BlanketOutboard 0, outer \n";
cout << (*BlanketOutboard).bounding_box_equatorial(1, "outer") << " BlanketOutboard 1, outer \n";
*/

if(blanketmode == "dcll") { // No OB HT Shield, placeholder = 0

 fileout << " Gap	;" 				<< (*BackWallOutboard).bounding_box_equatorial(0, "outer")-(*BlanketOutboard).bounding_box_equatorial(1, "outer") << "	;"
                                    << (*BackWallOutboard).bounding_box_equatorial(0, "outer") << "\n";
 fileout << " Back Wall OB	;" 		<< (*BackWallOutboard).bounding_box_equatorial(1, "outer")-(*BackWallOutboard).bounding_box_equatorial(0, "outer") << "	;"
                                    << (*BackWallOutboard).bounding_box_equatorial(1, "outer") << "\n";
 fileout << " Gap	;"				<< (*HTShieldDCLL).bounding_box_equatorial(1, "outer")-(*BackWallOutboard).bounding_box_equatorial(1, "outer") << "	;"
                                    << (*HTShieldDCLL).bounding_box_equatorial(1, "outer") << "\n";
 fileout << " HT Shield OB	;" 		<< (*HTShieldDCLL).bounding_box_equatorial(1, "outer")-(*HTShieldDCLL).bounding_box_equatorial(1, "inner") << "	;"
                                    << (*HTShieldDCLL).bounding_box_equatorial(1, "outer") << "\n"; 
 fileout << " Gap	;"				<< (*SkeletonRing).bounding_box_equatorial(1, "outer")-(*HTShieldDCLL).bounding_box_equatorial(1, "outer") << "	;"
									<< (*SkeletonRing).bounding_box_equatorial(1, "outer") << "\n";
 fileout << " Skeleton Ring OB	;"	<< (*SkeletonRing).bounding_box_equatorial(1, "outer")-(*SkeletonRing).bounding_box_equatorial(1, "inner") << "	;"
                                    << (*SkeletonRing).bounding_box_equatorial(1, "outer") << "\n";
 fileout << " Gap	;" 				<< (*VacuumVesselDCLL)[1].bounding_box_equatorial(0, "outer")-(*SkeletonRing).bounding_box_equatorial(1, "outer") << "	;"
                                    << (*VacuumVesselDCLL)[1].bounding_box_equatorial(0, "outer") << "\n";
 fileout << " Vacuum Vessel OB	;"	<< (*VacuumVesselDCLL)[1].bounding_box_equatorial(1, "outer")-(*VacuumVesselDCLL)[1].bounding_box_equatorial(0, "outer") << "	;"
                                    << (*VacuumVesselDCLL)[1].bounding_box_equatorial(1, "outer") << "\n";
 fileout << " Gap	;"				<< (*TFC).bounding_box_equatorial(1, "inner")-(*VacuumVesselDCLL)[1].bounding_box_equatorial(1, "outer") << "	;"
									<< (*TFC).bounding_box_equatorial(1, "inner") << "\n";
/*cout << (*BackWallOutboard).bounding_box_equatorial(0, "outer") << " BackWallOutboard 0, outer\n";
cout << (*BackWallOutboard).bounding_box_equatorial(1, "outer") << " BackWallOutboard 1, outer\n";
cout << (*HTShieldDCLL).bounding_box_equatorial(1, "outer") << " HTShieldDCLL 1, outer\n";
cout << (*HTShieldDCLL).bounding_box_equatorial(1, "inner") << " HTShieldDCLL 1, inner\n";
cout << (*SkeletonRing).bounding_box_equatorial(1, "outer") << " SkeletonRing 1, outer\n";
cout << (*SkeletonRing).bounding_box_equatorial(1, "inner") << " SkeletonRing 1, inner\n";
cout << (*VacuumVesselDCLL)[1].bounding_box_equatorial(0, "outer") << " (*VacuumVesselDCLL)[1] 0, outer\n";
cout << (*VacuumVesselDCLL)[1].bounding_box_equatorial(1, "outer") << " (*VacuumVesselDCLL)[1] 1, outer\n";
*/
}

else {
 fileout << " Gap	;" 				<< (*BlanketII).bounding_box_equatorial(0, "outer")-(*BlanketOutboard).bounding_box_equatorial(1, "outer") << "	;"
                                    << (*BlanketII).bounding_box_equatorial(0, "outer") << "\n";
 fileout << " Blanket II	;" 		<< (*BlanketII).bounding_box_equatorial(1, "outer")-(*BlanketII).bounding_box_equatorial(0, "outer") << "	;"
                                    << (*BlanketII).bounding_box_equatorial(1, "outer") << "\n";
 fileout << " Gap	;" 				<< (*HTShield)[1].bounding_box_equatorial(1, "inner")-(*BlanketII).bounding_box_equatorial(1, "outer") << "	;"
                                    << (*HTShield)[1].bounding_box_equatorial(1, "inner") << "\n";
 fileout << " HT Shield OB	;" 		<< (*HTShield)[1].bounding_box_equatorial(1, "outer")-(*HTShield)[1].bounding_box_equatorial(1, "inner") << "	;"
                                    << (*HTShield)[1].bounding_box_equatorial(1, "outer") << "\n";  // LCC HT Shield OK
 fileout << " Gap	;" 				<< (*VVessel)[1].bounding_box_equatorial(1, "inner")-(*HTShield)[1].bounding_box_equatorial(1, "outer") << "	;"
                                    << (*VVessel)[1].bounding_box_equatorial(1, "inner") << "\n";
 fileout << " Vacuum Vessel	;" 		<< (*VVessel)[1].bounding_box_equatorial(1, "outer")-(*VVessel)[1].bounding_box_equatorial(1, "inner") << "	;"
                                    << (*VVessel)[1].bounding_box_equatorial(1, "outer") << "\n";  // LCC VV OK
 fileout << " Gap	;" 				<< (*TFC).bounding_box_equatorial(1, "inner")-(*VVessel)[1].bounding_box_equatorial(1, "outer") << "	;"
                                    << (*TFC).bounding_box_equatorial(1, "inner") << "\n";
/*cout << (*BlanketII).bounding_box_equatorial(0, "outer") << " BlanketII 0, outer \n";
cout << (*BlanketII).bounding_box_equatorial(1, "outer") << " BlanketII 1, outer \n";
cout << (*HTShield)[1].bounding_box_equatorial(1, "outer") << " (*HTShield)[1] 1, outer \n";
cout << (*HTShield)[1].bounding_box_equatorial(1, "inner") << " (*HTShield)[1] 1, inner \n";
cout << (*VVessel)[1].bounding_box_equatorial(1, "outer") << " (*VVessel)[1] 1, outer \n";
cout << (*VVessel)[1].bounding_box_equatorial(1, "inner") << " (*VVessel)[1] 1, inner \n";
*/
}

// Combined
 fileout << " TF Coil	;" 			<< (*TFC).bounding_box_equatorial(1, "outer")-(*TFC).bounding_box_equatorial(1, "inner") << "	;"
                                    << (*TFC).bounding_box_equatorial(1, "outer") << "\n";  // LCC TF Coil OK
 fileout << " Gap	;" 				<< (*Cryo).bounding_box(2, "inner")-(*TFC).bounding_box_equatorial(1, "outer") << "	;"
                                    << (*Cryo).bounding_box(2, "inner") << "\n";
 fileout << " Cryostat	;"			<< (*Cryo).bounding_box(2, "outer")-(*Cryo).bounding_box(2, "inner") << "	;"
                                    << (*Cryo).bounding_box(2, "outer") << "\n";
/*cout << (*TFC).bounding_box_equatorial(1, "outer") << " TFC 1, outer \n";
cout << (*TFC).bounding_box_equatorial(1, "inner") << " TFC 1, inner \n";
cout << (*Cryo).bounding_box(2, "outer") << " Cryo 2, outer \n";
cout << (*Cryo).bounding_box(2, "inner") << " Cryo 2, inner \n";
*/

// ***** SA Testing *****
//cout << "Plasma SA using SA fcn    " << (*Plasma).surface_area("outer") <<"\n"; // 474 m2
//cout << "Plasma perimeter  " << (*Plasma).surface_area("plasma_perimeter") <<"\n"; // 14.3877 m
//cout << "Plasma SA using perim = " << (*Plasma).surface_area("plasma_perimeter")*2*pi*plasma_major_R <<"\n";
// special "plasma_perimeter" in Part.cpp surface_area calculation
// "*Part.perimeter" only works for Magnets in Cost.cpp, not for any other part
// ***** SA Testing *****


fileout << "\n\n";
fileout << "===============================================================\n";
fileout << "Geometric Data: Volumes [m^3]: \n\n";

// AT Plasma:
fileout << "Plasma: \n";
double Vplasma, Vaat, deltaV, Vplasma_SOL; // volumes of the systems code and ARIES-AT for comparison
if(blanketmode == "dcll") { 
		Vplasma = (*Plasma_DCLL).volume();
		fileout << " " << (*Plasma_DCLL).name() << ":	;" 	<< Vplasma << "\n\n"; }
else {	Vplasma = (*Plasma).volume();
		fileout << " " << (*Plasma).name() << ":	;" 		<< Vplasma << "\n\n"; }
// Vplasma_SOL = (*Plasma_AT_w_SOL).volume();
// Vaat = 308.219;  // AT plasma volume (95)
// fileout << " ARIES-AT Plasma:	;" 				<< Vaat << "\n";
// deltaV = 2.0*(Vplasma-Vaat)/(Vplasma+Vaat)*100.0;
// fileout << " Plasma Volume compared to AT:	;"	<< "[ " << deltaV << "% ]" << "\n\n";


// First Wall, LCC added 16x
fileout << "First Wall: \n";
 fileout << " " << (*FWInboard).name() << ":	;" 	<< 16.0*(*FWInboard).volume() << "\n";
 fileout << " " << (*FWOutboard).name() << ":	;" 	<< 16.0*(*FWOutboard).volume() << "\n";
 double Vfw = 16.0*((*FWInboard).volume() + (*FWOutboard).volume());
 fileout << "  Total First Wall Volume:	;" 			<< Vfw << "\n\n";
// Vaat = 3.234;  // ARIES-AT FW Vol, LCC removed /16
// deltaV = 2.0*(Vfw-Vaat)/(Vfw+Vaat)*100.0;
// fileout << "  ARIES-AT FW:	;" 						<< Vaat << "\n";
// fileout << "  FW Volume compared to AT:	;"			<< "[ " << deltaV << "% ]" << "\n\n";

 
// Blankets, LCC added 16x
// Aries-AT = 266.902 m3
fileout << "Blankets: \n";
 fileout << " " << (*BlanketInboard).name() << ":	;"	<< 16.0*(*BlanketInboard).volume() << "\n";
 fileout << " " << (*BlanketOutboard).name() << ":	;"	<< 16.0*(*BlanketOutboard).volume() << "\n";
double Vbl;
if(blanketmode == "dcll") {
 Vbl = 16.0*((*BlanketInboard).volume() + (*BlanketOutboard).volume());
 fileout << "  Total Blankets:	;"						<< Vbl << "\n\n"; }
else {
 Vbl = 16.0*((*BlanketInboard).volume() + (*BlanketOutboard).volume() + (*BlanketII).volume());
 fileout << " " << (*BlanketII).name() << ":	;"		<< 16.0*(*BlanketII).volume() << "\n";
 fileout << "  Total Blankets:	;"						<< Vbl << "\n\n"; }

// Back Wall (DCLL-only)
if(blanketmode == "dcll") {
fileout << "Back Wall: \n";
 fileout << " " << (*BackWallInboard).name() << ":	;"	<< 16.0*(*BackWallInboard).volume() << "\n";
 fileout << " " << (*BackWallOutboard).name() << ":	;"	<< 16.0*(*BackWallOutboard).volume() << "\n";
 fileout << "  Total Back Wall:	;"						<< 16.0*((*BackWallInboard).volume() + 
 														   		 (*BackWallOutboard).volume()) << "\n\n"; }


// HT Shields, LCC added 16x, ~176 m3, Xueren gets ~171
// 43. Penetration shielding fraction of HT shielding, read input file
  double pen_shielding, penetration_shielding;
  ifstream fin_pen;  fin_pen.open("InputData/powerflow.data");
  fin_pen.width(10); fin_pen.precision(2);
  fin_pen.setf(ios_base::uppercase, ios_base::scientific);
  int itxt_pen = 57; int inum_pen = 13; int iskip_pen = itxt_pen + inum_pen;
  pen_shielding; fin_pen.seekg(43*iskip_pen + itxt_pen +1); fin_pen >> pen_shielding;
	fin_pen.close();

fileout << "HT Shields: \n";
if(blanketmode == "dcll") { // no OB HT Shield
 fileout << " (No OB HT Shield for DCLL) \n";
 fileout << " " << (*HTShieldDCLL).name() << ":	;" 		<< 16.0*(*HTShieldDCLL).volume() << "\n";
 penetration_shielding = pen_shielding*16.0*(*HTShieldDCLL).volume();
 fileout << " Penetration Shielding Thickness:	;"	 	<< penetration_shielding << "\n";
 fileout << "  Total HT Shields:	;" 				  	<< 16.0*(*HTShieldDCLL).volume() + 
                                             			   penetration_shielding << "\n\n"; }
else {
 fileout << " " << (*HTShield)[0].name() << ":	;" 		<< 16.0*(*HTShield)[0].volume() << "\n";
 fileout << " " << (*HTShield)[1].name() << ":	;"		<< 16.0*(*HTShield)[1].volume() << "\n";
 penetration_shielding = pen_shielding*16.0*((*HTShield)[0].volume() + (*HTShield)[1].volume()); 
 fileout << " Penetration Shielding Thickness:	;"	 	<< penetration_shielding << "\n";
 fileout << "  Total HT Shields:	;" 				  	<< 16.0*((*HTShield)[0].volume() + 
                                             			   		 (*HTShield)[1].volume()) +
                                             			   penetration_shielding << "\n\n"; }

// Skeleton Ring (DCLL-only)
if(blanketmode == "dcll") {
fileout << "Skeleton Ring (IB+OB): \n";
 fileout << " " << (*SkeletonRing).name() 	<< ":	;" 	<< 16.0*(*SkeletonRing).volume() << "\n\n"; }


// Replaceable HT Shields, LCC added 16x 
fileout << "Replaceable HT Shields: \n";
fileout <<  " " << (*RepHTShield)[0].name() << ":	;" 	<< 16.0*( *RepHTShield)[0].volume() << "\n";
fileout <<  " " << (*RepHTShield)[1].name() << ":	;" 	<< 16.0*( *RepHTShield)[1].volume() << "\n";
fileout << "  Total Repl. HT Shields:	;"	 			<< 16.0*((*RepHTShield)[0].volume() + 
												  		   		 (*RepHTShield)[1].volume()) << "\n\n";


// Divertor Plates, LCC added 16x
fileout << "Divertor Plates (6): \n";
  double DivVol = 0.0;
  for(int i=0; i<6; i++)
  { fileout << " " << (*Divertor)[i].name() << ":	;"
                   << (*Divertor)[i].volume()*16.0 << "\n";
                   DivVol = DivVol + (*Divertor)[i].volume(); }
    fileout << "  Total Divertor:	;" <<  DivVol*16.0 << "\n\n"; // ~5.72 m3, Xueren gets 5.6 m3
            

// Vacuum Vessel, LCC added 16x, Aries-AT ~300 m3, Xueren ~393 m3, code ~267...
// This is the volume of the vessel itself, NOT the enclosed volume (est ~525 m3)
// Curiously "VacuumVesselMain.data" = "VacuumVesselOutboard.data"
// But break between main and OB occurs near midplane at think-to-thin interface
fileout << "Vacuum Vessel Volume (not enclosed volume): \n";
if(blanketmode == "dcll") {
fileout << " " << (*VacuumVesselDCLL)[0].name() << ":	;"	<< 16.0*( *VacuumVesselDCLL)[0].volume() << "\n";
fileout << " " << (*VacuumVesselDCLL)[1].name() << ":	;"	<< 16.0*( *VacuumVesselDCLL)[1].volume() << "\n";
fileout << " " << (*VacuumVesselDCLL)[2].name() << ":	;"	<< 16.0*( *VacuumVesselDCLL)[2].volume() << "\n";
fileout << " " << (*VacuumVesselDCLL)[3].name() << ":	;"	<< 16.0*( *VacuumVesselDCLL)[3].volume() << "\n";
fileout << "  Total Vacuum Vessel:	;" 						<< 16.0*((*VacuumVesselDCLL)[0].volume() + 
																	 (*VacuumVesselDCLL)[1].volume() +
															 		 (*VacuumVesselDCLL)[2].volume() +
															 		 (*VacuumVesselDCLL)[3].volume()) << "\n\n"; }
else {
fileout << " " << (*VVessel)[0].name() << ":	;"	<< 16.0*( *VVessel)[0].volume() << "\n";
fileout << " " << (*VVessel)[1].name() << ":	;"	<< 16.0*( *VVessel)[1].volume() << "\n";
fileout << "  Total Vacuum Vessel:	;" 				<< 16.0*((*VVessel)[0].volume() + 
															 (*VVessel)[1].volume()) << "\n\n"; }

// Maintenance Ports, should be 16 maintenance ports and vacuum ducts
// Duct 0,1,2 should comprise a volume = bottom+top+2*side
fileout << "Maintenance Ports and VV Ducts Volume (16) (not enclosed volume): \n";
fileout << " " << (*MPort).name() << ":	;"			<< 16.0*(*MPort).volume() << "\n";
fileout << " (MaintenancePort enclosed volume):	;"	<< 16.0*(*MPort).enclosed_volume() << "\n";
fileout << " " << (*Duct)[0].name() << ":	;"		<< 16.0*(*Duct)[0].volume() << "\n";
fileout << " " << (*Duct)[1].name() << ":	;"		<< 16.0*(*Duct)[1].volume() << "\n";
fileout << " " << (*Duct)[2].name() << "*2:	;"		<< 16.0*(*Duct)[2].volume()*2.0 << "\n";
double volMP = 16.0*((*MPort).volume()    +  (*Duct)[0].volume() + 
                     ( *Duct)[1].volume() + ((*Duct)[2].volume())*2.0);
fileout << "  Total All Ports and Ducts:	;"		<< volMP << "\n\n";
//cout << "( *MPort).volume() = " << ( *MPort).volume() << "\n";
            

// TF Coil, LCC added 16x
fileout << "TF Coils (16): \n";
fileout << " Total TF Coils:	;" 			<< 16.0*(*TFC).volume() << "\n\n";
 
 
// Other volumes: Bucking Cylinder, Caps, Cryo Dome, Manifolding, LCC added 16x
// LCC BCyl vol = 39.8 m^3, confirmed with CAD (i.e. need 16x)
// LCC CryoDome = 1041 m^3, comfirmed with CAD (i.e. need 16x)
fileout << "Other Volumes: \n";
fileout << " " << (*BCyl).name() << ":	;"	<< 16.0*(*BCyl).volume() << "\n";  
fileout << " TF Coil Upper Cap:	;" 			<< 16.0*(*UpperCap).volume() << "\n";
fileout << " TF Coil Lower Cap:	;" 			<< 16.0*(*LowerCap).volume() << "\n";
fileout << " Cryostat or Crydome:	;"	   	<< 16.0*(*Cryo).volume() << "\n";
fileout << " " << (*Manifolding)[0].name() << ":	;"	<< 16.0*(*Manifolding)[0].volume() << "\n";
fileout << " " << (*Manifolding)[1].name() << ":	;"	<< 16.0*(*Manifolding)[1].volume() << "\n\n";


// PF Coils including spares (both cost and volume)
// LCC added 16x to PF coil volumes since volume calculated by Geometry.cpp is UNIT vol
fileout << "PF Coils (36 active + 14 preliminary spares): \n";
double PoloFCVolume       = 0.0;
double PoloFCVolumeSpares = 0.0;
 for (int i=0; i<36; i++){ 
//   fileout << " " << (*PoloFC)[i].name() << ":	;" << (*PoloFC)[i].volume()*16.0 << "\n";
   PoloFCVolume = PoloFCVolume + (*PoloFC)[i].volume(); }  

// Spares included, #4-18
  for(int i=4; i<18; i++) 
  { PoloFCVolumeSpares = PoloFCVolumeSpares + (*PoloFC)[i].volume(); }
  fileout << "  Active PF Coil Vol. (36):	;"	<< 16.0*PoloFCVolume << "\n";  
  fileout << "  Spare  PF Coil Vol. (14):	;"	<< 16.0*PoloFCVolumeSpares << "\n";
  fileout << "   Total PF Coils Vol.:	;"		<< 16.0*(PoloFCVolume + PoloFCVolumeSpares) << "\n\n";
   
 
// Total First Wall, Blankets & Shielding (FWBS) and Fusion Power Core (FPC)
double FPC_vol = fpc_vol(); // calculated above, call DP.cpp ~780
if(blanketmode == "dcll") {
 fileout << " Total FW & Blankets:	;" <<           // include HT shielding???
                    ((*FWInboard).volume() +
                     (*FWOutboard).volume() +
                     (*BlanketInboard).volume() +
                     (*BlanketOutboard).volume() +
                     (*BackWallInboard).volume() +
                     (*BackWallOutboard).volume() +
                     (*SkeletonRing).volume())*16.0 << "\n";
 fileout << "  (IB & OB FW, Blanket & Back Wall, Skeleton Ring) \n\n"; 
 
 fileout << " Total Power Core (FPC) Vol:	;" << FPC_vol << "\n";
 fileout << "  (IB & OB FW, Blanket & Back Wall, Skeleton Ring, Divertor, \n"
            "   TF & PF magnets, VV, HT & Repl. Shield w/ Penetration Shielding) \n\n"; }
 
else {
 fileout << " Total FW & Blankets:	;" <<
                    ((*FWInboard).volume() +
                     (*FWOutboard).volume() +
                     (*BlanketInboard).volume() +
                     (*BlanketOutboard).volume() +
                     (*BlanketII).volume())*16.0 << "\n";
 fileout << "  (IB & OB FW & Blanket, Blanket II) \n\n";
 
 fileout << " Total Power Core (FPC) Vol:	;" << FPC_vol << "\n";
 fileout << "  (IB, OB FW & Blanket, Blanket II, Divertor, TF & PF magnets, \n"
            "   VV, HT & Repl. Shield w/ Penetration Shielding) \n\n"; }


// Support Structure Volume:
 fileout << " Support Structure Volume:	;" << 0.20*FPC_vol << "\n";
 fileout << "  (20% of FPC vol) \n\n"; // ~369 m3
//cout << "fpc_vol = " << FPC_vol << "\n";


// Reactor Building Volume:
 fileout << " Reactor Building Volume:	;" << reactor_building_volume() << "\n"; // from DP ~54
 fileout << "  (Uses exterior building dimensions) \n"; // ~140,000 m3
//cout << reactor_building_volume() << "\n"; 



fileout << "\n\n";
fileout << "===============================================================\n";
// Surface Areas
// LCC added 16x.  Hand calc w/ CAD, FW IB SA=136, FW OB SA=374
// ARIES-AT FW surface area = 425.588 m^2, ACT aggr = 474.711
// Multiple SA calc methods below:

double A_plasma_surface1, A_plasma_surface2, A_plasma_surface3;
if(blanketmode == "dcll")
		A_plasma_surface1 = (*Plasma_DCLL).surface_area("plasma_perimeter")*2.0*pi*plasma_major_R;
else 	A_plasma_surface1 = (*Plasma).surface_area("plasma_perimeter")*2.0*pi*plasma_major_R;

	A_plasma_surface2 = 4.0*pi*pi*plasma_major_R*(plasma_major_R/(*Outputs)[1])*
					    sqrt((1+pow((*Outputs)[12],2.0))/2.0); // elliptical method
	A_plasma_surface3 = 4.0*pi*pi*plasma_major_R*(plasma_major_R/(*Outputs)[1]+
					   	rad_gen_inb[0])*sqrt((1.0+pow((*Outputs)[12],2.0))/2.0);					   

 fileout << "Surface Areas [m^2]: \n\n";
if(blanketmode == "dcll") {
 fileout << " Plasma Surface Area:	;" 					<< (*Plasma_DCLL).surface_area("outer") << "\n"; 
 fileout << " Plasma Surface Area including SOL:	;" 	<< (*Plasma_AT_w_SOL).surface_area("outer") << "\n"; } // includes SOL, uses more accurate algorithm in G.cpp ~40
else { 
 fileout << " Plasma Surface Area:	;" 					<< (*Plasma).surface_area("outer") << "\n"; 
 fileout << " Plasma Surface Area including SOL:	;" 	<< (*Plasma_AT_w_SOL).surface_area("outer") << "\n"; } // includes SOL, uses more accurate algorithm in G.cpp ~40

 fileout << "  (using discretized geometry formulas)\n"; 
 fileout << " Scrape-off layer (SOL) thickness:	;"		<< rad_gen_inb[0] << "\n\n";

 //fileout << " Plasma Surface Area using perimeter method:	;" 	<< A_plasma_surface1 << "\n";
 //fileout << " Plasma Surface Area w/out SOL:	;" 		<< A_plasma_surface2 << "\n";
 //fileout << " Plasma Surface Area incl. SOL:	;" 		<< A_plasma_surface3 << "\n";
 
 
// FW Surface Area, why /2 ??
 fileout << " " << (*FWInboard).name() << ":	;"	<< 16.0/2.0*(*FWInboard).surface_area("outer") << "\n";
 fileout << " " << (*FWOutboard).name() << ":	;"	<< 16.0/2.0*(*FWOutboard).surface_area("outer") << "\n";
 double areaFW = 16.0/2.0*( (*FWInboard).surface_area("outer") + (*FWOutboard).surface_area("outer") );
 fileout << "  Total First Wall:	;" 				<< areaFW << "\n\n";

// Divertors: LCC added 16x
// ARIES-AT divertor surface area = 148.956 m^2, now ~191
 double DivArea = 0.0;
 for(int i=0; i<6; i++)
 { fileout << " " << (*Divertor)[i].name() << ":	;"
            	  << (*Divertor)[i].surface_area("outer")*16.0 << "\n";
            		 DivArea = DivArea + (*Divertor)[i].surface_area("outer")*16.0; }
 fileout << "  Total Divertors:	;" 		<<  DivArea << "\n";



fileout << "\n\n";
fileout << "===============================================================\n";
// Liquid Metal Mass
 fileout << "Liquid Metal Mass [tonne] \n\n";
 	// *16 for 1/16th segmentation, /1000 for kg to tonne
 double m_LM = 0.0;
 fileout << " " << (*FWInboard).name() << ":	;" 			<< (*FWInboard).liquid_metal_mass()*16.0/1000.0 << "\n";
 	m_LM = m_LM +  (*FWInboard).liquid_metal_mass(); // add up total LM mass
 fileout << " " << (*FWOutboard).name() << ":	;" 			<< (*FWOutboard).liquid_metal_mass()*16.0/1000.0 << "\n";
 	m_LM = m_LM +  (*FWOutboard).liquid_metal_mass();
 fileout << " " << (*BlanketInboard).name() << ":	;" 		<< (*BlanketInboard).liquid_metal_mass()*16.0/1000.0 << "\n";
 	m_LM = m_LM +  (*BlanketInboard).liquid_metal_mass();
 fileout << " " << (*BlanketOutboard).name() << ":	;" 		<< (*BlanketOutboard).liquid_metal_mass()*16.0/1000.0 << "\n";
 	m_LM = m_LM +  (*BlanketOutboard).liquid_metal_mass();

if(blanketmode == "dcll") {
 fileout << " " << (*BackWallInboard).name() << ":	;" 		<< (*BackWallInboard).liquid_metal_mass()*16.0/1000.0 << "\n";
 	m_LM = m_LM +  (*BackWallInboard).liquid_metal_mass();
 fileout << " " << (*BackWallOutboard).name() << ":	;" 		<< (*BackWallOutboard).liquid_metal_mass()*16.0/1000.0 << "\n";
 	m_LM = m_LM +  (*BackWallOutboard).liquid_metal_mass();
 fileout << " " << (*SkeletonRing).name() << ":	;" 			<< (*SkeletonRing).liquid_metal_mass()*16.0/1000.0 << "\n\n";
 	m_LM = m_LM +  (*SkeletonRing).liquid_metal_mass(); }
else {
 fileout << " " << (*BlanketII).name() << ":	;" 			<< (*BlanketII).liquid_metal_mass()*16.0/1000.0 << "\n\n";
 	m_LM = m_LM +  (*BlanketII).liquid_metal_mass(); }
 
 
// All divertor concepts now are He-cooled so LM_div is not needed
 double m_LM_div = 0.0;
/* for(int i=0; i<6; i++) 
 { fileout << " " << (*Divertor)[i].name() << ":	;"
            	  << (*Divertor)[i].liquid_metal_mass()*16.0/1000.0 << "\n";
   m_LM_div = m_LM_div + (*Divertor)[i].liquid_metal_mass(); }
   fileout << "  Total Divertors:	;" 				<< m_LM_div*16.0/1000.0 << "\n\n";
*/ 
 
double m_LM_rhts = 0.0;
 for(int i=0; i<2; i++) 
 { fileout << " " << (*RepHTShield)[i].name() << ":	;"
            	  << (*RepHTShield)[i].liquid_metal_mass()*16.0/1000.0 << "\n";
   m_LM_rhts = m_LM_rhts + (*RepHTShield)[i].liquid_metal_mass(); }
   fileout << "  Total Repl HT Shields:	;" 			<< m_LM_rhts*16.0/1000.0 << "\n\n";


double m_LM_man = 0.0;
 for(int i=0; i<2; i++) 
 { fileout << " " << (*Manifolding)[i].name() << ":	;"
            	  << (*Manifolding)[i].liquid_metal_mass()*16.0/1000.0 << "\n";
   m_LM_man = m_LM_man + (*Manifolding)[i].liquid_metal_mass(); }
   fileout << "  Total Manifolding:	;" 			<< m_LM_man*16.0/1000.0 << "\n\n";
 
 
double m_LM_hts = 0.0;
 if(blanketmode == "dcll") {
 fileout << " " << (*HTShieldDCLL).name() << ":	;" 	<< (*HTShieldDCLL).liquid_metal_mass()*16.0/1000.0 << "\n\n";
 	m_LM_hts = m_LM_hts + (*HTShieldDCLL).liquid_metal_mass(); }
 else {
 for(int i=0; i<2; i++) 
 { fileout << " " << (*HTShield)[i].name() << ":	;"
            	  << (*HTShield)[i].liquid_metal_mass()*16.0/1000.0 << "\n";
   m_LM_hts = m_LM_hts + (*HTShield)[i].liquid_metal_mass(); }
   fileout << "  Total HT Shields:	;" 				<< m_LM_hts*16.0/1000.0 << "\n\n"; }


// Total Liquid Metal Mass (AT was PbLi-cooled, 5269 tonnes total, 2108 in FWBS)
// "2.5 x power core vol = total coolant vol", L. Waganer, Total PbLi mass = mass of PbLi in FPC x 2.5
 m_LM = (m_LM + m_LM_div + m_LM_rhts + m_LM_man + m_LM_hts);
 fileout << "  Total Mass Liquid Coolant:	;" 		<< m_LM*16.0/1000.0 << "\n";
 fileout << "  Total w/ 2.5x BOP Factor:	;"		<< m_LM*16.0/1000.0*2.5 << "\n"; // 6312 tonnes



fileout << "\n\n";
fileout << "===============================================================\n";
// Helium Volume
 fileout << "Helium Volume [m^3] \n\n";
 fileout << "(Note that He density is taken as an avg between Brayton cycle temperatures and pressures.) \n";
if(blanketmode == "dcll") {
 fileout << " He vol for divertors, blanket and storage system, incl BOP: TBD	; \n"; }
else { 
 fileout << " He vol for divertors and storage system, incl BOP: TBD	; \n"; }

// include He in manifolding


fileout << "\n\n";
fileout << "===============================================================\n";
// Costs
 fileout << "Material-based Component Costs [M$]: \n";
 fileout << " See output file 'CostingAccounts.data' for the complete breakdown of all costing accounts.\n\n";
// LCC should these costs be taken from CA.cpp, inserted with account number???
// would have to call void CostingAccount::show() somehow

// First Walls and Blanket
 fileout << " " << (*FWInboard).name() << ":	;" 		<< (*FWInboard).cost()		 *16.0/1.e6 << "\n";
 fileout << " " << (*FWOutboard).name() << ":	;" 		<< (*FWOutboard).cost()		 *16.0/1.e6 << "\n";
 fileout << " " << (*BlanketInboard).name() << ":	;" 	<< (*BlanketInboard).cost()	 *16.0/1.e6 << "\n";
 fileout << " " << (*BlanketOutboard).name() << ":	;" 	<< (*BlanketOutboard).cost() *16.0/1.e6 << "\n";
double costfwbl;
if(blanketmode == "dcll") {
 fileout << " " << (*BackWallInboard).name() << ":	;" 	<< (*BackWallInboard).cost() *16.0/1.e6 << "\n";
 fileout << " " << (*BackWallOutboard).name() << ":	;" 	<< (*BackWallOutboard).cost()*16.0/1.e6 << "\n";
 fileout << " " << (*SkeletonRing).name() << ":	;" 		<< (*SkeletonRing).cost()    *16.0/1.e6 << "\n";
 costfwbl = ((*FWInboard).cost() 		 + (*FWOutboard).cost() +
             (*BlanketInboard).cost()  + (*BlanketOutboard).cost() +
             (*BackWallInboard).cost() + (*BackWallOutboard).cost() +
             (*SkeletonRing).cost())*16.0/1.e6;
 fileout << "  Total First Wall & Blankets:	;" << costfwbl << "\n\n"; }
else {
 fileout << " " << (*BlanketII).name() << ":	;" 		<< (*BlanketII).cost()		 *16.0/1.e6 << "\n"; 
 costfwbl = ((*FWInboard).cost() 		+ (*FWOutboard).cost() +
             (*BlanketInboard).cost() + (*BlanketOutboard).cost() +
             (*BlanketII).cost())*16.0/1.e6;
 fileout << "  Total First Wall & Blankets:	;" << costfwbl << "\n\n"; }
 
// Divertors
double DivCost = 0.0;
for(int i=0; i<6; i++) {
 fileout << " " << (*Divertor)[i].name() << ":	;" 		<< (*Divertor)[i].cost()*16.0/1.e6 << "\n";
 DivCost = DivCost + (*Divertor)[i].cost(); }
 fileout << "  Total Divertors (6):	;" 					<< DivCost*16.0/1.e6 << "\n\n";

 
// HT Shield
 fileout << " " << (*RepHTShield)[0].name() << ":	;" 	<< (*RepHTShield)[0].cost()	*16.0/1.e6 << "\n";
 fileout << " " << (*RepHTShield)[1].name() << ":	;" 	<< (*RepHTShield)[1].cost()	*16.0/1.e6 << "\n";
if(blanketmode == "dcll") {
 fileout << " " << (*HTShieldDCLL).name() << ":	;" 		<< (*HTShieldDCLL).cost()   *16.0/1.e6 << "\n";
 fileout << "  Total HT Shields:	;" <<   
 		  ((*RepHTShield)[0].cost() + (*RepHTShield)[1].cost() + (*HTShieldDCLL).cost())*16.0/1.e6 << "\n\n"; }
else {
 fileout << " " << (*HTShield)[0].name() << ":	;" 		<< (*HTShield)[0].cost()	*16.0/1.e6 << "\n";
 fileout << " " << (*HTShield)[1].name() << ":	;" 		<< (*HTShield)[1].cost()	*16.0/1.e6 << "\n"; 
 fileout << "  Total HT Shields:	;" <<   
 		  ((*RepHTShield)[0].cost() + (*RepHTShield)[1].cost() +
           (*HTShield)[0].cost() 	+ (*HTShield)[1].cost())*16.0/1.e6 << "\n\n"; }
 
// Vacuum Vessel
double VVcost = 0.0;
if(blanketmode == "dcll") {
for(int i=0; i<4; i++) {
 fileout << " " << (*VacuumVesselDCLL)[i].name() << ":	;" 	<< (*VacuumVesselDCLL)[i].cost()*16.0/1.e6 << "\n";
 VVcost = VVcost + (*VacuumVesselDCLL)[i].cost(); }
 fileout << "  Total Vacuum Vessel:	;" 					<< VVcost*16.0/1.e6 << "\n\n"; }
else {
 fileout << " " << (*VVessel)[0].name() << ":	;" 		<< (*VVessel)[0].cost()		*16.0/1.e6 << "\n";
 fileout << " " << (*VVessel)[1].name() << ":	;" 		<< (*VVessel)[1].cost()		*16.0/1.e6 << "\n";
 fileout << "  Total Vacuum Vessel:	;" <<  ((*VVessel)[0].cost() + (*VVessel)[1].cost())*16.0/1.e6 << "\n\n"; }
 
 
// Maintenance Ports, 16 maintenance ports and vacuum ducts
 fileout << " " << (*MPort).name() << ":	;"	<< 16.0*(*MPort).cost()/1.e6 << "\n";
 fileout << " " << (*Duct)[0].name() << ":	;"	<< 16.0*(*Duct)[0].cost()/1.e6 << "\n";
 fileout << " " << (*Duct)[1].name() << ":	;"	<< 16.0*(*Duct)[1].cost()/1.e6 << "\n";
 fileout << " " << (*Duct)[2].name() << "*2:	;"	<< 16.0*(*Duct)[2].cost()/1.e6 << "\n";
 double costMP = 16.0*((*MPort).cost()    + (*Duct)[0].cost() + 
          			   (*Duct)[1].cost() + ((*Duct)[2].cost()*2.0))/1.e6;
 fileout << "  Total Ports (16):	;"		 	<< costMP << "\n\n";
                                             			 
// TF Coil Cost
 fileout << " Total TF Coils Cost (16):	;" 		<< (*TFC).cost()*16/1.e6 << "\n\n"; // LCC 16x OK
 
 
// Bucking Cylinder, Caps, Cryo Dome, Manifolding
 fileout << " " << (*BCyl).name() << ":	;" 		<< (*BCyl).cost()	 *16.0/1.e6 << "\n";
 fileout << " TF Coil Upper Cap:	;" 			<< (*UpperCap).cost()*16.0/1.e6 << "\n";
 fileout << " TF Coil Lower Cap:	;" 			<< (*LowerCap).cost()*16.0/1.e6 << "\n";
 fileout << " " << (*Cryo).name() << ":	;" 		<< (*Cryo).cost()	 *16.0/1.e6 << "\n";
 fileout << " " << (*Manifolding)[0].name() << ":	;" << (*Manifolding)[0].cost()*16.0/1.e6 << "\n";
 fileout << " " << (*Manifolding)[1].name() << ":	;" << (*Manifolding)[1].cost()*16.0/1.e6 << "\n\n";


// PF Coil Cost, including spares
 fileout << " PF Coil Costs: 36 active, 14 spares (preliminary estimate):\n";
 double PoloFCCost       = 0.0;
 double PoloFCCostSpares = 0.0;
 for(int i=0; i<36; i++) { // was 29 but for what reason? artificially limited.
//  fileout << " " << (*PoloFC)[i].name() << ":	;" << (*PoloFC)[i].cost()*16.0/1.e6 << "\n";
  PoloFCCost = PoloFCCost + (*PoloFC)[i].cost(); }
// Spares included, called #4-18  
 for(int i=4; i<18; i++) {
  PoloFCCostSpares = PoloFCCostSpares + (*PoloFC)[i].cost(); }
  fileout << "  Active PF Coil Cost (36):	;"	<< PoloFCCost		*16.0/1.e6 << "\n";  
  fileout << "  Spare  PF Coil Cost (14):	;"	<< PoloFCCostSpares	*16.0/1.e6 << "\n";
  fileout << "   Total PF Coil Cost:	;"		<< PoloFCCost		*16.0/1.e6 +
  												   PoloFCCostSpares	*16.0/1.e6<< "\n";


fileout << "\n\n";
fileout << "===============================================================\n";
fileout << "Full Part Printouts: \n\n"; // new ~12/2011

if(blanketmode == "dcll") 	fileout << "-- DCLL blanket parts -- \n";
else						fileout << "-- SCLL blanket parts -- \n";

 fileout << "Part Name	;Material Name	;Volume Frac	;Density Frac	;Material Density	;Cost/Unit Mass	;Cost Basis Year	;Volume of Part	;Surface Area	;Liquid Metal Mass	;Mass of Part	;Cost of Part	;	\n";
 fileout << "	;	;	;	;[kg/m3]	;[$/kg]	; 	;[m3]	;[m2]	;[tonne]	;[kg]	;[M$]	;\n";
// ADD material description and purpose, also a summation row

/* This section was made into a separate function in Part.cpp ~996 for simplicity
// FWInboard
 int nmatcomp = (*FWInboard).number_material_components();
 fileout	<< (*FWInboard).name() << "	;";	
 for(int i=0; i<nmatcomp; i++) { // iterates # of materials, calls Part.cpp ~975
 fileout	<< (*FWInboard).part_components(i) << "	;";
  			if(i == nmatcomp-2)	fileout << "\n	;"; 
  			else { fileout  
 			<< 16.0*(*FWInboard).volume() << "	;"
 			<< 16.0/2.0*(*FWInboard).surface_area("outer") << "	;"
 			<< (*FWInboard).liquid_metal_mass()*16.0/1000.0 << "	;"
 			<< (*FWInboard).mass() << "	;"
 			<< (*FWInboard).cost()*16.0/1.e6 << "	;"; } }
 			fileout << "\n";
*/

// SCLL full part printouts (SCLL only for now...)
fileout	<< (*Plasma).name();	
fileout << "	;	;	;	;	;	;	;" << Vplasma << "	;" << (*Plasma_AT_w_SOL).surface_area("outer") << "\n"; // includes SOL

fileout	<< (*FWInboard).name();			fileout << (*FWInboard).part_printout();
fileout	<< (*FWOutboard).name();		fileout << (*FWOutboard).part_printout();
fileout << "  FW Sums	;	;	;	;	;	;	;" << Vfw << "	;" << areaFW << "\n";

fileout	<< (*BlanketInboard).name();	fileout << (*BlanketInboard).part_printout();
fileout	<< (*BlanketOutboard).name();	fileout << (*BlanketOutboard).part_printout();
fileout	<< (*BlanketII).name();			fileout << (*BlanketII).part_printout();
fileout << "  Blankets Sums	;	;	;	;	;	;	;" << Vbl << "\n";

fileout	<< (*HTShield)[0].name();		fileout << (*HTShield)[0].part_printout();
fileout	<< (*HTShield)[1].name();		fileout << (*HTShield)[1].part_printout();

fileout	<< (*RepHTShield)[0].name();	fileout << (*RepHTShield)[0].part_printout();
fileout	<< (*RepHTShield)[1].name();	fileout << (*RepHTShield)[1].part_printout();

fileout	<< (*Manifolding)[0].name();	fileout << (*Manifolding)[0].part_printout();
fileout	<< (*Manifolding)[1].name();	fileout << (*Manifolding)[1].part_printout();

for(int i=0; i<6; i++) {
fileout << (*Divertor)[i].name();		fileout << (*Divertor)[i].part_printout(); }
fileout << "  Divertors Sums	;	;	;	;	;	;	;" << DivVol*16.0 << "	;" << DivArea << "	;	;	;" << DivCost*16.0/1.e6  << "\n";

fileout	<< (*VVessel)[0].name();		fileout << (*VVessel)[0].part_printout();
fileout	<< (*VVessel)[1].name();		fileout << (*VVessel)[1].part_printout();

fileout	<< (*MPort).name();				fileout << (*MPort).part_printout();
fileout << "  MPort Sums	;	;	;	;	;	;	;" << volMP << "	;	;	;	;" << costMP  << "\n";

fileout	<< (*Duct)[0].name();			fileout << (*Duct)[0].part_printout();
fileout	<< (*Duct)[1].name();			fileout << (*Duct)[1].part_printout();
fileout	<< (*Duct)[2].name();			fileout << (*Duct)[2].part_printout();

fileout	<< (*Cryo).name();				fileout << (*Cryo).part_printout();
fileout	<< "TF Magnets (16)";			fileout << (*TFC).part_printout();
fileout	<< (*BCyl).name();				fileout << (*BCyl).part_printout();

fileout	<< "TF Coil Upper Cap";			fileout << (*UpperCap).part_printout();
fileout	<< "TF Coil Lower Cap";			fileout << (*LowerCap).part_printout();


for(int i=0; i<36; i++) {
fileout << (*PoloFC)[i].name();			fileout << (*PoloFC)[i].part_printout(); }
fileout << "  PF Coils Sums (#1-35)	;	;	;	;	;	;	;" << PoloFCVolume*16.0 << "	;	;	;	;" << PoloFCCost*16.0/1.e6  << "\n";
fileout << "  PF Coils Spares Sums (#4-18)	;	;	;	;	;	;	;" << PoloFCVolumeSpares*16.0 << "	;	;	;	;" << PoloFCCostSpares*16.0/1.e6  << "\n";


//double desc = (*FWInboard).material_description("SiC-comp");
//cout << desc << "\n";


fileout << "\n\n";
} // eof

// ====================================================================================
// END of detailed printout for "system" toggle, "void DesignPoint::show" function
// ====================================================================================


// ====================================================================================
// Extract Parameters from other modules

// Extract Plasma Parameters, from Aries.cpp ~200
  double DesignPoint::plasma_major_radius()			{ return plasma_params[0]; }
  double DesignPoint::plasma_minor_radius()			{ return plasma_params[1]; }
  double DesignPoint::plasma_elongation()			{ return plasma_params[2]; }
  double DesignPoint::plasma_triangularity()		{ return plasma_params[3]; }
  double DesignPoint::plasma_outer_squareness()		{ return plasma_params[4]; }
  double DesignPoint::plasma_inner_squareness()		{ return plasma_params[5]; }
  double DesignPoint::plasma_SOL()					{ return plasma_params[6]; }


// Extract the X-point of the last closed magnetic surface
  double DesignPoint::xpoint(int coordinate)
	{ if(coordinate == 0) return x_point_radius;
  	  if(coordinate == 1) return x_point_height;
  	  cout << "DesignPoint::xpoint(int coordinate): bad coordinate. \n";
      return 1.0e30; }


// TF coil inner radius of the inner leg and outer radius of the outer leg
  double DesignPoint::tf_magnet_inner_leg_radius() { return tf_coil_inboard_radius; }
  double DesignPoint::tf_magnet_outer_leg_radius() { return tf_coil_outboard_radius; }
  double DesignPoint::tf_magnet_max_field()		   { return tf_coil_max_magnetic_field; }
  int    DesignPoint::tf_magnet_cable_turns()      { return number_of_tf_cable_turns; }


// Extract Inboard Radial Builds, see "Read Radial Build" in Aries.cpp
	double DesignPoint::inboard_generic_build(int pick)
		{ for(int i=0; i<50; i++)
    	{ if(pick==i) {return rad_gen_inb[i]; } }
      	return 0.0; }

// Extract Outboard Radial Builds, see Aries.cpp, line ~228
	double DesignPoint::outboard_generic_build(int pick)
		{ for(int i=0; i<50; i++)
    	{ if(pick==i) {return rad_gen_outb[i];} }
     	return 0.0; }

// Extract Vertical Builds, see Aries.cpp, line ~230
	double DesignPoint::vertical_generic_build(int pick)
		{ for(int i=0; i<50; i++)
    	{ if(pick==i) {return top_gen_bld[i];} }
      	return 0.0; }

// Extract Neutron Wall Loads
double DesignPoint::neutron_wall_load(string nwl_id)
	{ if(nwl_id == "average" ) return ave_neutron_wall_load;
 	  if(nwl_id == "maximum" ) return max_neutron_wall_load;
	  cout << "DesignPoint::neutron_wall_load(string nwl_id): bad nwl_id. \n";
	  return 1.0e30; }


// Extract Maximum Heat Flux on Divertor
double DesignPoint::maximum_heat_flux_on_divertor(string location)
	{ if(location == "inboard")  return max_heat_flux_div_inboard;
	  if(location == "outboard") return max_heat_flux_div_outboard;
	  cout << "DesignPoint::maximum_heat_flux_on_divertor(string location): bad location. \n";
	  return 1.e30; }


// Extract Heat Fluxes on First Wall
double DesignPoint::fw_heat_flux(string flux_id)
	{ if(flux_id == "average" ) return ave_heat_flux_fw;
	  if(flux_id == "maximum" ) return max_heat_flux_fw;
	  cout << "DesignPoint::fw_heat_flux(string flux_id): bad flux_id. \n";
	  return 1.e30; }


// Extract Plasma Edge Safety Factor
double DesignPoint::plasma_edge_safety_factor() 
{ return q95; }

// Extract Plasma Current
double DesignPoint::plasma_current() 
{ return plsm_current; }

// Extract Plasma Poloidal Beta
double DesignPoint::plasma_poloidal_beta() 
{ return plsm_poloidal_beta; }

// FPC Volume, LCC added
double DesignPoint::fpc_vol() 
{ return fpcvolume; }

// Reactor building volume
double DesignPoint::reactor_building_volume() 
{ return vol_reactor_bldg; }

// Support Structure Volume
double DesignPoint::support_structure_volume() 
{ return vol_str; }


// ====================================================================================
// Current Drive 
// parameters are set in input file powerflow.data
// read from Aries.cpp ~532 "Aries_AT.get_cd_cost()"

void DesignPoint::get_cd_cost() {
//void DesignPoint::get_cd_cost(double Pow_NB, double Pow_LH, double Pow_FW, double freq_FW) {

// NEW - Read in CD parameters from powerflow.data input file
	ifstream fin_cd;  fin_cd.open("InputData/powerflow.data");
	fin_cd.width(10); fin_cd.precision(2);
	fin_cd.setf(ios_base::uppercase, ios_base::scientific);
	int itxt_cd = 57;  // length of first line of powerflow.data
	int inum_cd = 13;  // the number called is 13 letters long
	int iskip_cd = itxt_cd + inum_cd; // combination goes to the next parameter

// These are read in from input file but only FW is used.
	double Pow_NB;  fin_cd.seekg(36*iskip_cd + itxt_cd +1); fin_cd >> Pow_NB;
	double Pow_LH;  fin_cd.seekg(38*iskip_cd + itxt_cd +1); fin_cd >> Pow_LH;
	double Pow_FW;  fin_cd.seekg(40*iskip_cd + itxt_cd +1); fin_cd >> Pow_FW;
	double freq_FW; fin_cd.seekg(42*iskip_cd + itxt_cd +1); fin_cd >> freq_FW;

	double NB_efficiency; fin_cd.seekg(14*iskip_cd + itxt_cd +1); fin_cd >> NB_efficiency;
	double NB_unit_cost;  fin_cd.seekg(37*iskip_cd + itxt_cd +1); fin_cd >> NB_unit_cost;
	double LH_efficiency; fin_cd.seekg(15*iskip_cd + itxt_cd +1); fin_cd >> LH_efficiency;	
    double LH_unit_cost;  fin_cd.seekg(39*iskip_cd + itxt_cd +1); fin_cd >> LH_unit_cost;	
// EC???
	fin_cd.close();
	
// Set up table for ICRF fast wave (FW) efficiencies and costs
// Reference is old ASC Modules Docs "Power Flow.doc" by Raffray, Dragojlovic
// See also L. Waganer "Cost Basis Documentation (11-25-10).doc"
// Need to put these in an input file
    double frequency[6]  = {80.0e6, 158.0e6, 250.0e6, 800.0e6, 2500.0e6, 8000.0e6};
    double efficiency[6] = {0.84,   0.72,    0.65,    0.63,    0.48,     0.44    };
    double unit_cost[6]  = {1.65,   1.65,    3.29,    3.29,    3.29,     3.29    }; // 2009$
//    double unit_cost[6]  = {1.15,   1.15,    2.30,    2.30,    2.30,     2.87    }; // 1992$
    
    
// Interpolation from the CD table of frequencies for fast wave efficiency and unit cost
// Extrapolate frequencies WITHIN the defined domain
    for(int i=0; i<4; i++)
    { if(freq_FW >= frequency[i] && freq_FW < frequency[i+1])
		{ FW_efficiency = efficiency[i] +
                          (freq_FW - frequency[i]) *
                          (efficiency[i+1] - efficiency[i]) /
                          (frequency[i+1] - frequency[i]);
                                    
          FW_unit_cost =  unit_cost[i] +
                          (freq_FW - frequency[i]) *
                          (unit_cost[i+1] - unit_cost[i]) /
                          (frequency[i+1] - frequency[i]); } }
//cout << "freq FW (DP.cpp) = " << freq_FW << "\n"; // set to 95 MHz in input file
//cout << "CD efficiency (DP.cpp) = " << cd_fast_wave_efficiency << "\n"; // = 0.817


// Extrapolate frequencies BELOW the defined domain
    if(freq_FW < frequency[0])
       { FW_efficiency = efficiency[0] +
                         (freq_FW - frequency[0]) *
                         (efficiency[1] - efficiency[0]) /
                         (frequency[1] - frequency[0]);
                                   
          FW_unit_cost = unit_cost[0] +
                         (freq_FW - frequency[0]) *
                         (unit_cost[1] - unit_cost[0]) /
                         (frequency[1] - frequency[0]); }


// Extrapolate the frequencies ABOVE the defined domain
    if(freq_FW > frequency[5])
       { FW_efficiency = efficiency[4] +
                         (freq_FW - frequency[4]) *
                         (efficiency[5] - efficiency[4]) /
                         (frequency[5] - frequency[4]);
                                   
         FW_unit_cost = unit_cost[4] +
                        (freq_FW - frequency[4]) *
                        (unit_cost[5] - unit_cost[4]) /
                        (frequency[5] - frequency[4]); }


// CD costing
   cd_neutral_beam_total_cost = NB_unit_cost*Pow_NB;
   cd_lower_hybrid_total_cost = LH_unit_cost*Pow_LH;
   cd_fast_wave_total_cost    = FW_unit_cost*Pow_FW;
   cd_total_cost 			  = NB_unit_cost*Pow_NB + LH_unit_cost*Pow_LH + FW_unit_cost*Pow_FW;

} // eof
// ====================================================================================


// CD Costing
double DesignPoint::cd_cost()
{ return cd_total_cost; } // eof

// _ cost of the neutral beam plasma heating
double DesignPoint::cd_neutral_beam_cost()
{ return cd_neutral_beam_total_cost; } // eof

// _ total cost of the lower hybrid plasma heating
double DesignPoint::cd_lower_hybrid_cost()
{ return cd_lower_hybrid_total_cost; } // eof

// _ total cost of the fast wave plasma heating
double DesignPoint::cd_fast_wave_cost()
{ return cd_fast_wave_total_cost; } // eof


// Powers
double DesignPoint::gross_electrical_power()
{ return power_gross_electrical; } // eof

double DesignPoint::net_electrical_power()
{ return power_net_electrical; } // eof

double DesignPoint::thermal_power()
{ return power_thermal; } // eof

double DesignPoint::P_pump_blanket()
{ return power_blanket; } // eof

double DesignPoint::fusion_power()
{ return power_fusion; } // eof

double DesignPoint::data_point_filter(int loc) 
{ return data_flag[loc]; }
// ====================================================================================