// Geometry of Design Points 
// SiC blanket
// LCC cleaned up formatting ___ (still on-going)  Need to remove/archive old magnet algorithms

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

using namespace std;
#include <string>
using std::string;

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

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

// ====================================================================================
// 1/2 Plasma Geometry - at plasma surface
// Start with plasma, then build the other components around it.

// See ASC docs "Plasma Contour.doc" by Z.D.
// UPDATE: See  "Plasma Contour2 - Kessel - discretized formulas.doc"

void design_plasma_shape(DesignPoint *Machine, Part *Plasma)
	{ const double pi = 4.0*atan(1.0); // = 3.14159...    
      
// Plasma Parameters, from DP.cpp ~1540							// ACT-Ib
    const double Ro  = (*Machine).plasma_major_radius();		// 5.5
          double a   = (*Machine).plasma_minor_radius();		// 1.375
    const double k   = (*Machine).plasma_elongation();			// 2.2
    const double d   = (*Machine).plasma_triangularity();		// 0.7
    const double so  = (*Machine).plasma_outer_squareness(); 	// 1.0e-6
    const double si  = (*Machine).plasma_inner_squareness();	// 1.0e-6
    const double SOL = (*Machine).plasma_SOL();					// 0.1
    double theta     = 0.0;
    const double dt  = 0.01; // step size
    int N            = int(pi/dt+1);
    a = a; // add SOL to minor radius for plasma surface area calculation (see below)
      
// Evaluation of Contour Arrays
    int len = 2*N;
    double X[len], Y[len];
    double R, Z;

	for(int i=0; i<N; i++) // these are written to Aries_AT.data, ibl
		{ // ZD, uses asin(d)~d approx.
		  // R = Ro + a*cos(theta + d*sin(theta)); // not in use 
		  // CK, from Plasma Geometry Discretized Formulas.doc, more accurate
		    R = Ro + a*cos(theta + asin(d)*sin(theta)); 
	if(theta < pi/2.0)
        { Z = a*k*sin(theta + so*sin(2.0*theta)); }
	else
        { Z = a*k*sin(theta + si*sin(2.0*theta)); }
      	  X[N-i-1] = R;
	  	  Y[N-i-1] = Z;
      	  theta = theta + dt; // theta goes 0-pi by 0.01 
		//cout << theta << "	" << R << "	" << Z << "\n";
		}

// Close the loop by adding a symmetric profile in the clockwise direction:
	for(int i=0;  i<=N-1; i++)
    	{ X[N+i] =  X[N-i-1];
          Y[N+i] = -Y[N-i-1]; }


// Update the Plasma
	(*Plasma).get_contour_size(1,   "inner");
	(*Plasma).get_contour_size(len, "outer");

	double (*c_in)[2], (*c_out)[2], * d_in, * d_out;
	c_in  = new double [len][2];
	c_out = new double [len][2];
	d_in  = new double [len];
	d_out = new double [len];

	for(int i=0; i<len; i++)
	{ c_in[i][0]  = 0.0;
  	  c_in[i][1]  = 0.0;
  	  c_out[i][0] = X[i];
  	  c_out[i][1] = Y[i];
  	  d_in[i]     = 0.0;
  	  d_out[i]    = 2.0*pi*c_out[i][0]; }

	(*Plasma).get_contours(c_in, c_out, d_in, d_out);


// Volume and Surface Areas
	(*Plasma).get_volume();
	(*Plasma).get_surface_areas(); // calls Part.cpp ~641
//cout << "Plasma SA = " << (*Plasma).surface_area("outer") << "\n"; // ~472 m2 w/o SOL, ~503 w/ SOL

	delete [] c_in;
	delete [] c_out;
	delete [] d_in;
	delete [] d_out;
}   // eof


// ====================================================================================
// 2/2 Plasma Geometry - WITH SOL

// Copied from above, used for calculating plasma surface area w/ SOL in DP ~235
void design_plasma_shape_w_SOL(DesignPoint *Machine, Part *Plasma_AT_w_SOL) // Plasma_AT_w_SOL
	{ const double pi = 4.0*atan(1.0); // = 3.14159...    
      
// Plasma Parameters, from DP.cpp ~1540
    const double Ro  = (*Machine).plasma_major_radius();
          double a   = (*Machine).plasma_minor_radius();
    const double k   = (*Machine).plasma_elongation();
    const double d   = (*Machine).plasma_triangularity();
    const double so  = (*Machine).plasma_outer_squareness();
    const double si  = (*Machine).plasma_inner_squareness();
    const double SOL = (*Machine).plasma_SOL();	
    double theta     = 0.0;
    const double dt  = 0.01; // step size
    int N            = int(pi/dt+1);
    a = a + SOL; // ADD SOL to minor radius for plasma surface area calculation

// Evaluation of Contour Arrays
    int len = 2*N;
    double X[len], Y[len];
	double R, Z;

	for(int i=0; i<N; i++) // these are written to Aries_AT.data, ibl
		{ //R = Ro + a*cos(theta + d*sin(theta)); 	  // ZD, uses asin(d)~d approx.
		  // CK, from Plasma Geometry Discretized Formulas.doc, more accurate
		  R = Ro + a*cos(theta + asin(d)*sin(theta));
	if(theta < pi/2.0)
        { Z = a*k*sin(theta + so*sin(2.0*theta)); }
	else
        { Z = a*k*sin(theta + si*sin(2.0*theta)); }
      	  X[N-i-1] = R;
	  	  Y[N-i-1] = Z;
      	  theta = theta + dt; // theta goes 0-pi by 0.01 
		//cout << theta << "	" << R << "	" << Z << "\n";
		}

// Close the loop by adding a symmetric profile in the clockwise direction:
	for(int i=0;  i<=N-1; i++)
    	{ X[N+i] =  X[N-i-1];
          Y[N+i] = -Y[N-i-1]; }


// Update the Plasma
	(*Plasma_AT_w_SOL).get_contour_size(1,   "inner");
	(*Plasma_AT_w_SOL).get_contour_size(len, "outer");

	double (*c_in)[2], (*c_out)[2], * d_in, * d_out;
	c_in  = new double [len][2];
	c_out = new double [len][2];
	d_in  = new double [len];
	d_out = new double [len];

	for(int i=0; i<len; i++)
	{ c_in[i][0]  = 0.0;
  	  c_in[i][1]  = 0.0;
  	  c_out[i][0] = X[i];
  	  c_out[i][1] = Y[i];
  	  d_in[i]     = 0.0;
  	  d_out[i]    = 2.0*pi*c_out[i][0]; }

	(*Plasma_AT_w_SOL).get_contours(c_in, c_out, d_in, d_out);


// Volume and Surface Areas
	(*Plasma_AT_w_SOL).get_volume();
	(*Plasma_AT_w_SOL).get_surface_areas(); // calls Part.cpp ~641
//cout << "Plasma SA w/ SOL = " << (*Plasma_AT_w_SOL).surface_area("outer") << "\n"; // ~472 m2 w/o SOL, ~503 w/ SOL

	delete [] c_in;
	delete [] c_out;
	delete [] d_in;
	delete [] d_out;
}   // eof


// ====================================================================================
// Last Closed Magnetic Surface (LCMSurf or Separatrix)

void design_separatrix_shape(DesignPoint *Machine, Part * LCMSurf, Part * Plasma)
{       const double pi = 4.0*atan(1.0);
        double delta = (*Machine).inboard_generic_build(0)/2.0; // distance between plasma surface and separatrix = 1/2 SOL
        double dist  = 2.0;

// The angles alphain and alphaout are taken from CAD drawings of ARIES-AT                                                                
        double alphain  =      66.25*pi/180.0; // = 1.156 rad
	    double alphaout = pi - 19.44*pi/180;   // = 2.802 rad

// Plasma bounding box:
        double ksilo = (*Plasma).bounding_box(0, "outer");
        double ksihi = (*Plasma).bounding_box(2, "outer");

// Reference points:
        double xs[10], ys[10];
        
/* ... The X-point is determined as an intersection between the straight portions
   ... of the last closed magnetic surface on the inner and the outer side. */
       double xleft, yleft, kleft, xright, yright, kright;
        
// ... first reference point ...
        xs[0] = ksilo - delta;
        ys[0] = 0.0;


/* ... The second point is at the distance delta from the plasma contour, at the
	   location where the tangent to the contour forms a given angle with the r-axis.

   ... The third point is located on the straight line that extends from the second
       point at the angle alphain. */

      int lplasm =  (*Plasma).contour_size("outer");
      for(int i=0; i<lplasm; i++)
      {
      double xphi = (*Plasma).contour(i+1, 0, "outer");
      double yphi = (*Plasma).contour(i+1, 1, "outer");
      
      double xpmid = (*Plasma).contour(i, 0, "outer");
      double ypmid = (*Plasma).contour(i, 1, "outer");
      
      double xplo = (*Plasma).contour(i-1, 0, "outer");
      double yplo = (*Plasma).contour(i-1, 1, "outer");

	  double tagm = sqrt(pow((xphi-xplo),2) + pow((yphi-yplo),2));
      double alpha;
        
      if(ypmid > 0.0 && yphi > yplo)
        { alpha = asin((yphi - yplo)/tagm);
          if(alpha >= alphain)
          { double nx = -(yphi - yplo)/tagm;
	        double ny =  (xphi - xplo)/tagm;
      		double tx =  ny;
      		double ty = -nx;
      		
       xs[1] = xpmid + nx*delta;
	   ys[1] = ypmid + ny*delta;
       xs[2] = xpmid + tx*dist;
	   ys[2] = ypmid + ty*dist;
	  
// ... starting point and slop of the straight section
       xleft = xs[1]; yleft = ys[1]; kleft = tan(alphain); } }
} // end loop around plasma contour


     double ref_pts_i[10][2];
     ref_pts_i[0][0] =  xs[2];
     ref_pts_i[0][1] = -ys[2];
     
     ref_pts_i[1][0] =  xs[1];
     ref_pts_i[1][1] = -ys[1];
     
     ref_pts_i[2][0] =  xs[0];
     ref_pts_i[2][1] =  ys[0] - 0.01;
     
     ref_pts_i[3][0] =  xs[0];
     ref_pts_i[3][1] =  ys[0] + 0.01;
     
     ref_pts_i[4][0] =  xs[1];
     ref_pts_i[4][1] =  ys[1];
     
     ref_pts_i[5][0] =  xs[2];
     ref_pts_i[5][1] =  ys[2];

     int nptsi = 6;

      
// ... ... resolution ... ...
    int resol_i[100];
    for(int i=0; i<nptsi; i++)
    { resol_i[i] = 50; }
      resol_i[nptsi-1] = 0;


// ... ... connectors ... ...
    string connectors_i[100];
    connectors_i[0] = "line";
    connectors_i[1] = "target";
    connectors_i[2] = "line";
    connectors_i[3] = "target";
    connectors_i[4] = "line";
    connectors_i[5] = "none";
    

// ... Outer Loop ... reference points ... ...
     double ref_pts_o[100][2];

      xs[2] = ksihi + delta;
      ys[2] = 0.0;

      for(int i=0; i<lplasm; i++)
      { double xphi  = (*Plasma).contour(i+1, 0, "outer");
        double yphi  = (*Plasma).contour(i+1, 1, "outer");

        double xpmid = (*Plasma).contour(i, 0, "outer");
        double ypmid = (*Plasma).contour(i, 1, "outer");

        double xplo  = (*Plasma).contour(i-1, 0, "outer");
        double yplo  = (*Plasma).contour(i-1, 1, "outer");

        double tagm  = sqrt( pow( (xphi-xplo),2) + pow( (yphi-yplo),2) );
        double alpha;

      if( ypmid > 0.0 )
      { alpha = asin( (yphi-yplo)/tagm );
        if( alpha >= (alphaout-pi) )
        { double nx = -(yphi - yplo)/tagm;
	  	  double ny =  (xphi - xplo)/tagm;
      	  double tx =  ny; 
	  	  double ty = -nx;

      xs[1] = xpmid + nx*delta;
  	  ys[1] = ypmid + ny*delta;

      xs[0] = xs[1] - tx*dist;
  	  ys[0] = ys[1] - ty*dist;


// ... starting point and slope of the straight section
       xright = xs[1]; yright = ys[1]; kright = tan(alphaout); } }
}   // end loop around plasma contour


         ref_pts_o[0][0] =  xs[0];
         ref_pts_o[0][1] =  ys[0];
         
         ref_pts_o[1][0] =  xs[1];
         ref_pts_o[1][1] =  ys[1];
         
         ref_pts_o[2][0] =  xs[2];
         ref_pts_o[2][1] =  ys[2] - 0.01;
         
         ref_pts_o[3][0] =  xs[2];
         ref_pts_o[3][1] =  ys[2] + 0.01;
              
         ref_pts_o[4][0] =  xs[1];
         ref_pts_o[4][1] = -ys[1];
         
         ref_pts_o[5][0] =  xs[0];
         ref_pts_o[5][1] = -ys[0];

         int nptso = 6;


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }
      resol_o[nptso-1] = 0;


// ... ... connectors ... ...
    string connectors_o[100];
    connectors_o[0] = "line";
    connectors_o[1] = "target";
    connectors_o[2] = "line";
    connectors_o[3] = "target";
    connectors_o[4] = "line";
    connectors_o[5] = "none";


// ... Calculate the location of the X-point above the plasma:
     double x_xp = (yleft - yright + kright*xright - kleft*xleft)/(kright - kleft);
     double y_xp =  yleft + kleft*x_xp - kleft*xleft;
     (*Machine).get_xpoint(x_xp, y_xp);
    

// ... Update the last closed magnetic surface ...
(*LCMSurf).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i,
									ref_pts_o, connectors_i, connectors_o);

// Define the Inner Target Curve For Blanket Shape
	double targcurve[10000][2];
    for(int i=0; i<lplasm; i++)
    { targcurve[i][0] = (*Plasma).contour(i, 0, "outer");
      targcurve[i][1] = (*Plasma).contour(i, 1, "outer"); }


// Connect the Reference Points into Contours
 double center_i[100][2], center_o[100][2];
 (*LCMSurf).connect(center_i, center_o, targcurve, lplasm);


// Establish Contour Depth for volume calculation
 double * d_in; double * d_out;
 int lin  = (*LCMSurf).contour_size("inner");
 int lout = (*LCMSurf).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 2.0*pi*(*LCMSurf).contour(i, 0, "inner")/16.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*LCMSurf).contour(i, 0, "outer")/16.0;

 (*LCMSurf).get_depth(d_in, d_out);


// Volume and Surface Areas
(*LCMSurf).get_volume();
(*LCMSurf).get_surface_areas();

delete [] d_in;
delete [] d_out;
} // eof


// ====================================================================================
// First Wall IB Geometry

void design_first_wall_inboard_shape(DesignPoint * Machine, Part * Wall, Part * Plasma)
{ 
if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour ...
/* ... ... start from the outer bounding box of plasma as a reference entity and build on top of it ... */
     double ksilo = (*Plasma).bounding_box(0, "outer"); 
     double etahi = (*Plasma).bounding_box(3, "outer"); // height of plasma


// ... ... reference points ... ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];

     double dscr  = (*Machine).inboard_generic_build(0); // SOL thickness
     double dwall = (*Machine).inboard_generic_build(1); // FW thickness
     double scale = 0.877; // constant from First Wall.doc to scale FW height geometrically to AT


// Define the Outer Contour reference points ... ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 4;
     double xs[10], ys[10];

     xs[0] = ksilo - dscr - dwall;
     ys[0] = 0.0;

     xs[1] = xs[0];
     ys[1] = scale*etahi; 

     xs[2] = ksilo - dscr;
     ys[2] = ys[1];

     xs[3] = xs[2];
     ys[3] = 0.0;

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

/* ... ... Close the loop by adding a symmetric profile in the clockwise direction: ... ... */
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-2;


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }


// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line"; }


// LCC DCLL differs here, calls "blanket" the "first wall" instead but the equations are the same...

// Update the Blanket
(*Wall).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                 connectors_i, connectors_o);


// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2]; int tgl=0;


// Connect the Reference Points into Contours
(*Wall).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*Wall).contour_size("outer");
 
 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)  d_in[i]  = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Wall).contour(i, 0, "outer")/16.0;

 (*Wall).get_depth(d_in, d_out);


// Volume and Surface Areas
(*Wall).get_volume();
(*Wall).get_surface_areas();

delete [] d_in;
delete [] d_out;
} // eof


// ====================================================================================
// First Wall OB Geometry - SCLL (DCLL code is different but build looks ok - should implement DCLL)

void design_first_wall_outboard_shape(DesignPoint * Machine, Part * Wall, Part * Plasma)
{									// Ref Part is plasma
if( (*Machine).data_point_filter(0) > 0 ) return; 

// ... Define the Inner Contour ...
// start from outer bounding box of plasma as ref entity and build on top of it
     double ksilo = (*Plasma).bounding_box(0, "outer"); // IB plasma location
     double ksihi = (*Plasma).bounding_box(2, "outer"); // OB plasma location
     double etahi = (*Plasma).bounding_box(3, "outer"); // height of plasma
//cout << "etahi plasma = " << etahi << "\n";
     

// ... ... reference points ... ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];
     
     double dscr  = (*Machine).inboard_generic_build(0); // SOL thickness
     double dwall = (*Machine).inboard_generic_build(1); // FW thickness
	 double scale = 0.34833; // from ASC Docs First Wall.doc to scale geometrically to AT - 
	 // the point at which the FW turns upward to meet the repl HT shield 

// LCC: DCLL differs here....................

// Define the Outer Contour
// ... reference points ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 6;
     double xs[10], ys[10];

	// 1 Start of contour at midplane
     xs[0] = ksihi + dscr;
     ys[0] = 0.0;

	// 2 Point at which FW turns upward from following plasma
     xs[1] = ksilo - dscr + scale*(ksihi-ksilo); 
     ys[1] = 0.99*etahi; // nearly the height of the plasma, ref?

	// 3 Top of FW
     xs[2] = xs[1];
     ys[2] = etahi + 0; // + 1.2577; // ref? dist above plasma to extend FW to repl HT shield
	 // top of FW = top of plasma + 1.2577 = bottom of repl HT (used below in repl HT)
	
	// 4 Top-right of FW
     xs[3] = xs[2] + dwall;
     ys[3] = ys[2];

	// 5 Inside inflection point
     xs[4] = xs[3];
     ys[4] = ys[1] + dwall;

	// 6 End at midplane
     xs[5] = xs[0] + dwall;
     ys[5] = 0.0;

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

/* ... ... Close the loop by adding a symmetric profile in the clockwise direction: ... ... */
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0]; 
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
       nptso = 2*nptso-2;


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }


// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
      if(i == 0 || i == 4 || i == 5 || i == 9) connectors_o[i] = "target"; }


// Update the Blanket
(*Wall).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                 connectors_i, connectors_o);


// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2];
int tgl = (*Plasma).contour_size("outer");
    for(int i=0; i<tgl; i++)
    { targcurve[i][0] = (*Plasma).contour(i, 0, "outer");
      targcurve[i][1] = (*Plasma).contour(i, 1, "outer"); }


// Connect the Reference Points into Contours
(*Wall).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*Wall).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)  d_in[i]  = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Wall).contour(i, 0, "outer")/16.0;

 (*Wall).get_depth(d_in, d_out);


// Volume and Surface Areas
(*Wall).get_volume();
(*Wall).get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Blanket IB Geometry

void design_blanket_inboard_shape(DesignPoint * Machine, Part * Blanket, Part * RefBody)
{ 								// RefBody is FWInboard
if( (*Machine).data_point_filter(0) > 0 ) return; 

// ... Define the Inner Contour ...
// ... start from the first wall as the reference entity and build on top of it
     double ksilo = (*RefBody).bounding_box(0, "outer");
     double etahi = (*RefBody).bounding_box(3, "outer");

// ... reference points ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];
     
     double gap   = 0.0; 								 // gap between the FW and blanket
     double dwall = (*Machine).inboard_generic_build(2); // blanket thickness

// Define the Outer Contour
// ... reference points ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 4;
     double xs[10], ys[10];

     xs[0] = ksilo - gap - dwall; // removed xs[0] = ksilo - (1.0 - 0.04)*dwall;
     ys[0] = 0.0;

     xs[1] = xs[0];
     ys[1] = etahi;

     xs[2] = ksilo - gap; // removed xs[2] = ksilo;
     ys[2] = ys[1];

     xs[3] = xs[2];
     ys[3] = 0.0;

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

// ... Close the loop by adding a symmetric profile in the clockwise direction: ...
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
       nptso = 2*nptso-2;


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }


// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line"; }


// Update the Blanket
(*Blanket).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                 	connectors_i, connectors_o);

// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2]; int tgl=0;


// Connect the Reference Points into Contours
(*Blanket).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*Blanket).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)  d_in[i]  = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Blanket).contour(i, 0, "outer")/16.0;

 (*Blanket).get_depth(d_in, d_out);

// Volume and Surface Areas
(*Blanket).get_volume();
(*Blanket).get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Back Wall IB Geometry (DCLL)

void design_back_wall_inboard_shape(DesignPoint * Machine, Part * BackWall, Part * RefBody)
{ 									// RefBody is BlanketInboard
if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour ...
// ... start from the first wall as the reference entity and build on top of it
     double ksilo = (*RefBody).bounding_box(0, "outer");
     double etahi = (*RefBody).bounding_box(3, "outer");

// ... reference points ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];
     
     double gap   = (*Machine).inboard_generic_build(5); // blanket and back wall gap
     double dwall = (*Machine).inboard_generic_build(6); // back wall thickness


// Define the Outer Contour
// ... reference points ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 4;
     double xs[10], ys[10];

     xs[0] = ksilo - gap - dwall;
     ys[0] = 0.0;

     xs[1] = xs[0];
     ys[1] = etahi;

     xs[2] = ksilo -gap;
     ys[2] = ys[1];

     xs[3] = xs[2];
     ys[3] = 0.0;

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

//... Close the loop by adding a symmetric profile in the clockwise direction: ...
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
       nptso = 2*nptso-2;

// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line"; }


// Update the Back Wall
(*BackWall).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                     connectors_i, connectors_o);

// Define the Inner Target Curve For Back Wall Shape
double targcurve[10000][2]; int tgl=0;


// Connect the Reference Points into Contours
(*BackWall).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*BackWall).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*BackWall).contour(i, 0, "outer")/16.0;

 (*BackWall).get_depth(d_in, d_out);


// Volume and Surface Areas
(*BackWall).get_volume();
(*BackWall).get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Back Wall OB Geometry (DCLL)

void design_back_wall_outboard_shape(DesignPoint * Machine, Part * BackWall, Part * RefBody)
{									// RefBody is BlanketOutboard 
if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour ...
//... start from the first wall as the reference entity and build on top of it
     double ksilo = (*RefBody).bounding_box(0, "outer");
     double ksihi = (*RefBody).bounding_box(2, "outer");
     double etahi = (*RefBody).bounding_box(3, "outer");

// ... reference points ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];
     double dblanket = (*Machine).outboard_generic_build(2); // blanket I thickness
     double gap      = 0.0;									 // blanket I to back wall II gap
     double dwall    = (*Machine).outboard_generic_build(9); // back wall II thickness


// Define the Outer Contour
// ... reference points ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 4;
     double xs[10], ys[10];

      xs[0] = ksihi + gap;
      ys[0] = 0.0;
      
      xs[1] = ksilo + dblanket + gap;
      ys[1] = etahi;
      
      xs[2] = xs[1] + dwall;
      ys[2] = etahi;
      
      xs[3] = xs[0] + dwall;
      ys[3] = 0.0;

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

// ... Close the loop by adding a symmetric profile in the clockwise direction: 
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] = ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-2;

// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
      if(i == 0 || i == 2 || i == 3 || i == 5) connectors_o[i] = "target"; }


// Update the Blanket
(*BackWall).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                     connectors_i, connectors_o);

// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2];
int tgl = (*RefBody).contour_size("outer");
    for(int i=0; i<tgl; i++)
    { targcurve[i][0] = (*RefBody).contour(i, 0, "outer");
      targcurve[i][1] = (*RefBody).contour(i, 1, "outer"); }


// Connect the Reference Points into Contours
(*BackWall).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*BackWall).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*BackWall).contour(i, 0, "outer")/16.0;

 (*BackWall).get_depth(d_in, d_out);


// Volume and Surface Areas
(*BackWall).get_volume();
(*BackWall).get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Blanket OB Geometry

void design_blanket_outboard_shape(DesignPoint * Machine, Part * Blanket, Part * RefBody)
{ 									// RefBody is FWOutboard
if( (*Machine).data_point_filter(0) > 0 ) return; 

// ... Define the Inner Contour ...
// start from the first wall as the reference entity and build on top of it
     double ksilo = (*RefBody).bounding_box(0, "outer");
     double ksihi = (*RefBody).bounding_box(2, "outer");
     double etahi = (*RefBody).bounding_box(3, "outer"); // FWOutboard height
//cout << "etahi FW = " << etahi << "\n";

// ... ... reference points ... ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];
     
     double dfw   = (*Machine).outboard_generic_build(1); // FW I thickness
     double gap   = 0.0;								  // FW I to blanket gap
	 double dwall = (*Machine).outboard_generic_build(2); // OB blanket I thickness
	 // removed dwall = (*Machine).inboard_radial_build(3)*0.85714 (scaling????)

// Define the Outer Contour
// ... ... reference points ... ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 4;
     double xs[10], ys[10];

	// 1 OB midplane start
      xs[0] = ksihi + gap; // added +gap
      ys[0] = 0.0;

	// 2 Top of plasma where divertor meets      
      xs[1] = ksilo + dfw + gap; // removed xs[1] = ksilo + 0.0467*dwall
      ys[1] = etahi; // -0.1; // arbitrary

	// 3 Top of blanket where it meets repl HT shield and HT shield
      xs[2] = xs[1] + dwall; // removed xs[2] = xs[1] + (1.0-0.0467)*dwall
      ys[2] = etahi + 0.15; // -0.7; // arbitrary, sets angle

	// 4 Follow contour back down at thickness dwall to midplane      
      xs[3] = xs[0] + dwall; // removed xs[3] = xs[0] + (1.0-0.0467)*dwall
      ys[3] = 0.0;


     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

/* ... ... Close the loop by adding a symmetric profile in the clockwise direction: ... ... */
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-2;


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }


// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
      if(i == 0 || i == 2 || i == 3 || i == 5) connectors_o[i] = "target"; }


// Update the Blanket
(*Blanket).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                 	connectors_i, connectors_o);


// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2];
int tgl = (*RefBody).contour_size("outer");

    for(int i=0; i<tgl; i++)
    { targcurve[i][0] = (*RefBody).contour(i, 0, "outer");
      targcurve[i][1] = (*RefBody).contour(i, 1, "outer"); }


// Connect the Reference Points into Contours
(*Blanket).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*Blanket).contour_size("outer");
 
 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Blanket).contour(i, 0, "outer")/16.0;

 (*Blanket).get_depth(d_in, d_out);


// Volume and Surface Areas
(*Blanket).get_volume();
(*Blanket).get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Blanket II (SCLL only)

void design_blanket_II_shape(DesignPoint * Machine, Part * Blanket, Part * RefBody)
{ 											// RefBody = BlanketOutboard
if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour ...
// ... start from the FW as the reference entity and build on top of it ... 	
	 double ksilo = (*RefBody).bounding_box(0, "outer"); // inner OB blanket
     double ksihi = (*RefBody).bounding_box(2, "outer"); // outer OB blanket
     double etahi = (*RefBody).bounding_box(3, "outer"); // top of OB blanket
//cout << "ksilo = " << ksilo << "\n";
//cout << "ksihi = " << ksihi << "\n";
//cout << "etahi = " << etahi << "\n";

// ... ... reference points ... ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];
     double dwref = (*Machine).outboard_generic_build(2); // OB blkt thickness, rest goes to FW to sum = 35 cm
     double dwgap = (*Machine).outboard_generic_build(6); // OB blkt and blkt II gap (stabilizing shell?)
     double dwall = (*Machine).outboard_generic_build(8); // blkt II thickness
// Removed *0.85714 from (3) and *1.875 from (5), ALSO changed IB to OB radial build
//cout << "ob thick = " << dwref << "\n";
//cout << "ob II thick = " << dwall << "\n";

// Define the Outer Contour
// ... ... reference points ... ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 4;
     double xs[10], ys[10];

	// OB midplane start
      xs[0] = ksihi + dwgap;
      ys[0] = 0.0;

      xs[1] = ksilo + dwref + dwgap;
      ys[1] = etahi; // mod 1/31/12 for manifolding area

      xs[2] = xs[1] + dwall;
      ys[2] = etahi + 0.215; // to make angle from OB blnkt // sets angle, mod 1/31/12 for manifolding area

      xs[3] = xs[0] + dwall;
      ys[3] = 0.0;

/*cout << "xs[0], ys[0] = " << xs[0] << ",	" << ys[0] << "\n";
cout << "xs[1], ys[1] = " << xs[1] << ",	" << ys[1] << "\n";
cout << "xs[2], ys[2] = " << xs[2] << ",	" << ys[2] << "\n";
cout << "xs[3], ys[3] = " << xs[3] << ",	" << ys[3] << "\n";
*/

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

/* ... ... Close the loop by adding a symmetric profile in the clockwise direction: ... ... */
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-2;


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }


// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
      if(i == 0 || i == 2 || i == 3 || i == 5) connectors_o[i] = "target"; }


// Update the Blanket
(*Blanket).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                 	connectors_i, connectors_o);


// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2];
int tgl = (*RefBody).contour_size("outer");
    for(int i=0; i<tgl; i++)
    { targcurve[i][0] = (*RefBody).contour(i, 0, "outer");
      targcurve[i][1] = (*RefBody).contour(i, 1, "outer"); }


// Connect the Reference Points into Contours
(*Blanket).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*Blanket).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Blanket).contour(i, 0, "outer")/16.0;

 (*Blanket).get_depth(d_in, d_out);


// Volume and Surface Areas
(*Blanket).get_volume();
(*Blanket).get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Divertor Plates - DCLL (slightly different than SiC, especially IB side)
// Only difference was LHS ref was 1/3*

void design_divertor_shape_dcll(DesignPoint * Machine, Part (*Divertor)[6], Part * LCMSurf,
                           		Part * RefLeft, Part * RefRight, Part * Plasma)
						// LCMSurf = sepatrix, RefLeft = IB back wall, RefRight = OB FW
{
if( (*Machine).data_point_filter(0) > 0 ) return;

// Bounding Box
     double ksiloLeft 	= (*RefLeft).bounding_box(0, "outer");  // IB back wall
     double ksihiLeft 	= (*RefLeft).bounding_box(2, "outer");
     double etahiLeft 	= (*RefLeft).bounding_box(3, "outer");
     double ksiloRight 	= (*RefRight).bounding_box(0, "outer"); // OB FW
     double etahiRight 	= (*RefRight).bounding_box(3, "outer");
     double etahiPlsm 	= (*Plasma).bounding_box(3, "outer");
     
     double delta 		= (*Machine).inboard_generic_build(0)/10.0; // 1/10th SOL
     double dfw 		= (*Machine).inboard_generic_build(1); // FW thickness
     double dgfwb 		= 0.0;								   // FW to blanket gap
     double dblkt 		= (*Machine).inboard_generic_build(2); // blanket thickness
     double dgbbw 		= (*Machine).inboard_generic_build(5); // blanket to back wall gap thickness
     double dbw 		= (*Machine).inboard_generic_build(6); // back wall thickness
     double drephts 	= (*Machine).vertical_generic_build(3);// repl HT shield thickness

// Implement the curvature identical to the ARIES-AT divertor plates by stretching
// the normalized vectors which sample this curvature between the calculated
// end-points of the divertor plates generated by this code.

// Normalized sample vectors: (reference?)
double xnorm0[10] = {0.0,  0.01,  0.01,  0.03, -0.02, -0.17, -0.35, -0.56, -0.76, -1.0};
double ynorm0[10] = {0.0,  0.17,  0.31,  0.45,  0.57,  0.69,  0.79,  0.86,  0.93,  1.0};

double xnorm1[10] = {0.0,  0.09,  0.24,  0.37,  0.51,  0.67,  0.81,  0.92,  0.98,  1.0};
double ynorm1[10] = {0.0,  0.00,  0.00,  0.03,  0.10,  0.23,  0.39,  0.56,  0.75,  1.0};

double xnorm2[10] = {0.0,  0.0,   0.0,   0.0,   0.0,   0.0,   0.01,  0.19,  0.49,  1.0};
double ynorm2[10] = {0.0, -0.09, -0.18, -0.31, -0.43, -0.55, -0.66, -0.79, -0.90, -1.0};

double Xtransl, Ytransl, Cx, Cy; // transformation parameters


// Loop Over  Divertor Plates
for(int k=0; k<6; k++)
{ double ref_pts_o[100][2], center_o[100][2];
  int nptso = 10; // number of reference points
  double xs[10], ys[10];
  switch(k) // switch to determine reference points
{

// Reference Points for Plate 1
case 0 :
      xs[0] = ksihiLeft + dgbbw + dblkt + dgfwb + dfw; // top of IB FW
      ys[0] = etahiLeft;
      xs[9] = ksiloLeft + drephts;

for(int i=0; i<(*LCMSurf).contour_size("outer"); i++)
{ if( (*LCMSurf).contour(i, 0, "outer") > xs[9] &&
      (*LCMSurf).contour(i, 1, "outer") > 0.0 )
      { ys[9] = (*LCMSurf).contour(i, 1, "outer");
        break; } }


// Transformation from normalized to actual divertor plate:
Xtransl = xs[0]; Ytransl = ys[0];
for(int i=1; i<9; i++)
	{ xs[i] = xnorm0[i]*(xs[0]-xs[9]) + Xtransl;
      ys[i] = ynorm0[i]*(ys[9]-ys[0]) + Ytransl; } 
      break; // from case 0
      
      
// Reference points for Plate 2
case 1:
    xs[0] =  ksiloLeft + drephts;
    for(int i=0; i<(*LCMSurf).contour_size("outer"); i++)
{ if( (*LCMSurf).contour(i, 0, "outer") > xs[0] )
   	{ ys[0] = (*LCMSurf).contour(i, 1, "outer") + delta;
      break; } }
    ys[9] = etahiRight - drephts;
        for(int i=0; i<(*LCMSurf).contour_size("inner"); i++)
{ if( (*LCMSurf).contour(i, 1, "inner") < ys[9] )
    { xs[9] = (*LCMSurf).contour(i, 0, "inner"); } }


// Transformation from normalized to actual divertor plate:
Xtransl = xs[0]; Ytransl = ys[0];
for(int i=1; i<9; i++)
	{ xs[i] = xnorm1[i]*(xs[9]-xs[0]) + Xtransl;
 	  ys[i] = ynorm1[i]*(ys[9]-ys[0]) + Ytransl; }    
    break; // from case 1

  
// Reference points for Plate 3
case 2 :    
      ys[0] = etahiRight - drephts;
        for(int i=0; i<(*LCMSurf).contour_size("inner"); i++)
{ if( (*LCMSurf).contour(i, 1, "inner") < ys[0] )
	{ xs[0] = (*LCMSurf).contour(i, 0, "inner") + delta; } }
    xs[9] = (*RefRight).reference_point(1, 0, "outer");
    ys[9] = (*RefRight).reference_point(1, 1, "outer");


// Transformation from normalized to actual divertor plate:
Xtransl = xs[0]; Ytransl = ys[0];
for(int i=1; i<9; i++)
{ xs[i] = xnorm2[i]*(xs[9]-xs[0]) + Xtransl;
  ys[i] = ynorm2[i]*(ys[0]-ys[9]) + Ytransl; }  
    break; // from case 2


// Reference points for Plate 4  (they mirror Plate 3)
case 3 :
	xs[0] =  (*RefRight).reference_point(1, 0, "outer");
    ys[0] = -(*RefRight).reference_point(1, 1, "outer");
    ys[9] = etahiRight - drephts;
        
for(int i=0; i<(*LCMSurf).contour_size("inner"); i++)
{ if( (*LCMSurf).contour(i, 1, "inner") < ys[9] )
	{ xs[9] = (*LCMSurf).contour(i, 0, "inner") + delta; } }
      ys[9] = -ys[9];

// Transformation from normalized to actual divertor plate:
Xtransl = xs[0]; Ytransl = ys[0];
for(int i=1; i<9; i++)
	{ xs[i] = (xnorm2[9-i]-1.0)*(xs[0]-xs[9]) + Xtransl;
	  ys[i] = (-ynorm2[9-i]-1.0)*(ys[0]-ys[9]) + Ytransl; }    
     break; // from case 3


// Reference points for Plate 5  (they mirror Plate 2)
case 4 :
     ys[0] = etahiRight - drephts;
        for(int i=0; i<(*LCMSurf).contour_size("inner"); i++)
{ if( (*LCMSurf).contour(i, 1, "inner") < ys[0] )
	{ xs[0] = (*LCMSurf).contour(i, 0, "inner"); } }
      ys[0] = - ys[0];
      xs[9] = ksiloLeft + drephts;
    for(int i=0; i<(*LCMSurf).contour_size("outer"); i++)
{ if( (*LCMSurf).contour(i, 0, "outer") > xs[9] )
	{ ys[9] = (*LCMSurf).contour(i, 1, "outer") + delta;
      break; } }
      ys[9] = -ys[9];

// Transformation from normalized to actual divertor plate:
Xtransl = xs[9]; Ytransl = ys[9];
for(int i=1; i<9; i++)
	{ xs[i] =  xnorm1[9-i]*(xs[0]-xs[9]) + Xtransl;
      ys[i] = -ynorm1[9-i]*(ys[9]-ys[0]) + Ytransl; } 
      break; // from case 4 


// Reference points for Plate 6  (they mirror Plate 1)
case 5 :
      xs[0] = ksiloLeft + drephts;
	for(int i=0; i<(*LCMSurf).contour_size("outer"); i++)
{ if( (*LCMSurf).contour(i, 0, "outer") > xs[0] &&
      (*LCMSurf).contour(i, 1, "outer") > 0.0 )
       { ys[0] = -(*LCMSurf).contour(i, 1, "outer");
         break; } }
      xs[9] = ksihiLeft + dgbbw + dblkt + dgfwb + dfw;
      ys[9] = -etahiLeft;

// Transformation from normalized to actual divertor plate:
Xtransl = xs[9]; Ytransl = ys[9];
for(int i=1; i<9; i++)
	{ xs[i] =  xnorm0[9-i]*(xs[9]-xs[0]) + Xtransl;
      ys[i] = -ynorm0[9-i]*(ys[9]-ys[0]) + Ytransl; }  
      break; // from case 5
} // end switch over reference points

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

// ... ... resolution ... ...
    int resol_o[100];
     resol_o[0] = 10;
     resol_o[1] = 10;
     resol_o[2] = 10;
     resol_o[3] = 10;
     resol_o[4] = 10;
     resol_o[5] = 10;
     resol_o[6] = 10;
     resol_o[7] = 10;
     resol_o[8] = 10;
     resol_o[9] = 0;

// ... ... connectors ... ...
    string connectors_o[100];
    connectors_o[0] = "line";
    connectors_o[1] = "line";
    connectors_o[2] = "line";
    connectors_o[3] = "line";
    connectors_o[4] = "line";
    connectors_o[5] = "line";
    connectors_o[6] = "line";
    connectors_o[7] = "line";
    connectors_o[8] = "line";
    connectors_o[9] = "none";

// Update the Divertor Plate
int nptsi = 1;
int resol_i[100];
double ref_pts_i[100][2];
string connectors_i[100];

(*Divertor)[k].get_contour_descriptions(nptsi, nptso, 	resol_i, 	resol_o,
                                 		ref_pts_i, 		ref_pts_o,
                                 		connectors_i, 	connectors_o);

// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2];
int tgl = 0;

// ... Connect the Reference Points into Contours ...
double center_i[100][2];
(*Divertor)[k].connect(center_i, center_o, targcurve, tgl);

// ... Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin  = (*Divertor)[k].contour_size("inner");
 int lout = (*Divertor)[k].contour_size("outer");
 d_in  = new double [lin];
 d_out = new double [lout];
 for(int i=0; i<lin; i++)  d_in[i]  = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Divertor)[k].contour(i, 0, "outer")/16.0;
 (*Divertor)[k].get_depth(d_in, d_out);


// ... Volume and Surface Areas ...
double ddiv = (*Machine).vertical_generic_build(1); // divertor plate thickness
	(*Divertor)[k].get_surface_areas();
	(*Divertor)[k].get_shell_volume(ddiv); // shell volume multiplies the volume by thickness

delete [] d_in;
delete [] d_out;

} // end loop over divertor plates
} // eof


// ====================================================================================
// Divertor Plates - SCLL, there are 6 of these: upper and lower IB, OB, dome

void design_divertor_shape(DesignPoint * Machine, Part (*Divertor)[6], Part * LCMSurf,
                           Part * RefLeft, Part * RefRight, Part * Plasma)
						// LCMSurf = sepatrix, RefLeft = IB blanket, RefRight = OB blanket

{ if( (*Machine).data_point_filter(0) > 0) return;

// Bounding Box
     double ksiloLeft  = (*RefLeft).bounding_box( 0, "outer"); // IB blanket
     double ksihiLeft  = (*RefLeft).bounding_box( 2, "outer");
     double etahiLeft  = (*RefLeft).bounding_box( 3, "outer");
     double ksiloRight = (*RefRight).bounding_box(0, "outer"); // OB blanket
     double etahiRight = (*RefRight).bounding_box(3, "outer");
     double etahiPlsm  = (*Plasma).bounding_box(  3, "outer"); // top of plasma
     
     double delta      = (*Machine).inboard_generic_build(0)/10.0; // 1/10th SOL
     double dwall      = (*Machine).inboard_generic_build(2);      // IB blanket thickness
     double drephts    = (*Machine).vertical_generic_build(3);     // repl HT shield thickness
     double newtopref  = etahiPlsm + 1.2577 - drephts; // to replace OB blanket ref (etahi2)
      // Need to ref top of plasma + 1.2577 factor taken from OB FW part description
      // since divertors can no longer ref OB blanket due to manifolding, mod 1/30/12
      // top of divertors = top of plasma + 1.2577 - repl HT thickness = bottom of repl HT
      
// Implement the curvature identical to the ARIES-AT divertor plates by stretching the
// normalized vectors which sample this curvature between the calculated end-points
// of the divertor plates generated by this code.

// Normalized sample vectors: ref? ZD says divertors are represented by "Bezier curves"
// to allow the user to prescribe an arbitrary angle between div and LCMS...
double xnorm0[10] = {0.0,  0.01,  0.01,  0.03, -0.02, -0.17, -0.35, -0.56, -0.76, -1.0};
double ynorm0[10] = {0.0,  0.17,  0.31,  0.45,  0.57,  0.69,  0.79,  0.86,  0.93,  1.0};

double xnorm1[10] = {0.0,  0.09,  0.24,  0.37,  0.51,  0.67,  0.81,  0.92,  0.98,  1.0};
double ynorm1[10] = {0.0,  0.00,  0.00,  0.03,  0.10,  0.23,  0.39,  0.56,  0.75,  1.0};

double xnorm2[10] = {0.0,  0.0,   0.0,   0.0,   0.0,   0.0,   0.01,  0.19,  0.49,  1.0};
double ynorm2[10] = {0.0, -0.09, -0.18, -0.31, -0.43, -0.55, -0.66, -0.79, -0.90, -1.0};

double Xtransl, Ytransl, Cx, Cy; // transformation parameters


// Loop Over 6 Divertor Plates ----------------------
for(int k=0; k<6; k++)
{ double ref_pts_o[100][2], center_o[100][2];
  int nptso = 10; // number of reference points
  double xs[10], ys[10];
  switch(k) // switch to determine reference points
{
// Reference Points for Plate 1
    case 0 :
      xs[0] = ksihiLeft; // edge of IB blanket
      ys[0] = etahiLeft; // top of IB blanket
      xs[9] = ksiloLeft + drephts; // removed + (1.0/3.0)*dwall; // end at repl HT shield wall
      // 1/3*dwall was an estimate for repl HT thickness, now in radial build input file

for(int i=0; i<(*LCMSurf).contour_size("outer"); i++)
{  if( (*LCMSurf).contour(i, 0, "outer") > xs[9] &&
       (*LCMSurf).contour(i, 1, "outer") > 0.0 )
       { ys[9] = (*LCMSurf).contour(i, 1, "outer");
         break; } }


// Transformation from normalized to actual divertor plate:
Xtransl = xs[0]; Ytransl = ys[0];
for(int i=1; i<9; i++)
{   xs[i] = xnorm0[i]*(xs[0]-xs[9]) + Xtransl;
    ys[i] = ynorm0[i]*(ys[9]-ys[0]) + Ytransl; } 
    break; // from case 0
      
      
// Reference points for Plate 2
    case 1: 
    xs[0] =  ksiloLeft + drephts; // removed + (1.0/3.0)*dwall; // start of plate 2
    for(int i=0; i<(*LCMSurf).contour_size("outer"); i++)
{   if( (*LCMSurf).contour(i, 0, "outer") > xs[0] )
       { ys[0] = (*LCMSurf).contour(i, 1, "outer") + delta;
         break; } }

    ys[9] = newtopref; // etahiRight - drephts; // removed - 0.5*dwall; // end of plate 2
        for(int i=0; i<(*LCMSurf).contour_size("inner"); i++)
{   if( (*LCMSurf).contour(i, 1, "inner") < ys[9] )
       { xs[9] = (*LCMSurf).contour(i, 0, "inner"); } }


// Transformation from normalized to actual divertor plate:
Xtransl = xs[0]; Ytransl = ys[0];
for(int i=1; i<9; i++)
	{ xs[i] = xnorm1[i]*(xs[9]-xs[0]) + Xtransl;
	  ys[i] = ynorm1[i]*(ys[9]-ys[0]) + Ytransl; }   
    break; // from case 1
    

// Reference points for Plate 3
    case 2 :
      ys[0] = newtopref; // etahiRight - drephts; // removed - 0.5*dwall; // start of plate 3
        for(int i=0; i<(*LCMSurf).contour_size("inner"); i++)
{   if( (*LCMSurf).contour(i, 1, "inner") < ys[0] )
       { xs[0] = (*LCMSurf).contour(i, 0, "inner") + delta; } }
    xs[9] = ksiloRight;
    ys[9] = 0.99*etahiPlsm;


// Transformation from normalized to actual divertor plate:
Xtransl = xs[0]; Ytransl = ys[0];
for(int i=1; i<9; i++)
	{ xs[i] = xnorm2[i]*(xs[9]-xs[0]) + Xtransl;
      ys[i] = ynorm2[i]*(ys[0]-ys[9]) + Ytransl; } 
    break; // from case 2
    

// Reference points for Plate 4  (they mirror Plate 3)
    case 3 :
       xs[0] = ksiloRight;
       ys[0] = -0.99*etahiPlsm;
       ys[9] = newtopref; // etahiRight - drephts; // removed - 0.5*dwall;

for(int i=0; i<(*LCMSurf).contour_size("inner"); i++)
{   if( (*LCMSurf).contour(i, 1, "inner") < ys[9] )
       { xs[9] = (*LCMSurf).contour(i, 0, "inner") + delta; } }
       ys[9] = -ys[9];


// Transformation from normalized to actual divertor plate:
Xtransl = xs[0]; Ytransl = ys[0];
for(int i=1; i<9; i++)
	{ xs[i] =  (xnorm2[9-i]-1.0)*(xs[0]-xs[9]) + Xtransl;
  	  ys[i] = (-ynorm2[9-i]-1.0)*(ys[0]-ys[9]) + Ytransl; } 
     break; // from case 3


// Reference points for Plate 5  (they mirror Plate 2)
    case 4 :
     ys[0] = newtopref; // etahiRight - drephts; // removed - 0.5*dwall;
        for(int i=0; i<(*LCMSurf).contour_size("inner"); i++)
{   if( (*LCMSurf).contour(i, 1, "inner") < ys[0] )
       { xs[0] = (*LCMSurf).contour(i, 0, "inner"); } }
    ys[0] = - ys[0];
    xs[9] = ksiloLeft + drephts; // removed + (1.0/3.0)*dwall;
    for(int i=0; i<(*LCMSurf).contour_size("outer"); i++)
{   if( (*LCMSurf).contour(i, 0, "outer") > xs[9] )
       { ys[9] = (*LCMSurf).contour(i, 1, "outer") + delta;
         break; } }
    ys[9] = -ys[9];


// Transformation from normalized to actual divertor plate:
Xtransl = xs[9]; Ytransl = ys[9];
for(int i=1; i<9; i++)
	{ xs[i] =  xnorm1[9-i]*(xs[0]-xs[9]) + Xtransl;
      ys[i] = -ynorm1[9-i]*(ys[9]-ys[0]) + Ytransl; } 
    break; // from case 4


// Reference points for Plate 6  (they mirror Plate 1)
    case 5 :
      xs[0] = ksiloLeft + drephts; // removed + (1.0/3.0)*dwall;
for(int i=0; i<(*LCMSurf).contour_size("outer"); i++)
{   if( (*LCMSurf).contour(i, 0, "outer") > xs[0] &&
        (*LCMSurf).contour(i, 1, "outer") > 0.0 )
       { ys[0] = -(*LCMSurf).contour(i, 1, "outer");
         break; } }
      xs[9] = ksihiLeft;
      ys[9] = -etahiLeft;


// Transformation from normalized to actual divertor plate:
Xtransl = xs[9]; Ytransl = ys[9];
for(int i=1; i<9; i++)
	{ xs[i] =  xnorm0[9-i]*(xs[9]-xs[0]) + Xtransl;
      ys[i] = -ynorm0[9-i]*(ys[9]-ys[0]) + Ytransl; } 
    break; // from case 5


} // end switch over reference points
     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }


// ... ... resolution ... ...
    int resol_o[100];
     resol_o[0] = 10;     resol_o[1] = 10;
     resol_o[2] = 10;     resol_o[3] = 10;
     resol_o[4] = 10;     resol_o[5] = 10;
     resol_o[6] = 10;     resol_o[7] = 10;
     resol_o[8] = 10;     resol_o[9] = 0;
      

// ... ... connectors ... ...
    string connectors_o[100];
    connectors_o[0] = "line";    connectors_o[1] = "line";
    connectors_o[2] = "line";    connectors_o[3] = "line";
    connectors_o[4] = "line";    connectors_o[5] = "line";
    connectors_o[6] = "line";    connectors_o[7] = "line";
    connectors_o[8] = "line";    connectors_o[9] = "none";


// Update the Divertor Plate
	int nptsi = 1;
	int resol_i[100];
	double ref_pts_i[100][2];
	string connectors_i[100];

(*Divertor)[k].get_contour_descriptions(nptsi, 			nptso, 		resol_i,
										resol_o,		ref_pts_i, 	ref_pts_o,
		                                connectors_i, 	connectors_o);


// Define the Inner Target Curve For Blanket Shape
	double targcurve[10000][2];
	int tgl = 0;


// ... Connect the Reference Points into Contours ...
	double center_i[100][2];
	(*Divertor)[k].connect(center_i, center_o, targcurve, tgl);


// ... Establish Contour Depth for volume calculation ...
	const double pi = 4.0*atan(1.0); // pi = 3.14159....
 	double * d_in; 
 	double * d_out;
 	int lin  = (*Divertor)[k].contour_size("inner");
 	int lout = (*Divertor)[k].contour_size("outer");
 
 	d_in  = new double [lin];
 	d_out = new double [lout];
 
 	for(int i=0; i<lin; i++)  d_in[i]  = 0.0;
	for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Divertor)[k].contour(i, 0, "outer")/16.0;
	(*Divertor)[k].get_depth(d_in, d_out);


// ... Volume and Surface Areas ...
	double ddiv = (*Machine).vertical_generic_build(1); // divertor plate thickness, 3.35 cm (AT), ~7 cm (ACT)
	(*Divertor)[k].get_surface_areas();
	(*Divertor)[k].get_shell_volume(ddiv); // shell volume multiplies the volume by thickness
//cout << "(*Divertor) = " << (*Divertor)[0].volume() << "\n"; // assume vol returned is shell vol but hard to check

	delete [] d_in;
	delete [] d_out;

} // end loop over divertor plates
} // eof


// ====================================================================================
// Replaceable HT Shield, 2 parts: top and bottom (both blankets)

void design_replaceable_ht_shield(DesignPoint *Machine, Part (*Shield)[2],	Part *RefLo, 
										Part *RefHi, 	Part *Plasma)
							// RefLo = BackWallInboard   (DCLL), BlanketInboard  (SiC)
							// RefHi = FirstWallOutboard (DCLL), BlanketOutboard (SiC)
{ if( (*Machine).data_point_filter(0) > 0 ) return;

// Define the Inner Contour of Top Shield...
// ... start from the first wall as the reference entity and build on top of it
     double ksilo1 	  = (*RefLo).bounding_box(0, "outer");  // inside IB blanket
     double etahi1 	  = (*RefLo).bounding_box(3, "outer");  // top IB blanket
     double ksilo2 	  = (*RefHi).bounding_box(0, "outer");  // inside OB blanket
     double etahi2 	  = (*RefHi).bounding_box(3, "outer");  // top OB blanket
     double etahiplsm = (*Plasma).bounding_box(3, "outer"); // top plasma
      
/*cout << "ksilo1  = " << ksilo1 << "\n";
cout << "etahi1  = " << etahi1 << "\n";
cout << "ksilo2  = " << ksilo2 << "\n";
cout << "etahi2  = " << etahi2 << "\n";
cout << "etahiplsm  = " << etahiplsm << "\n";
*/

// ... ... reference points ... ...
    double ref_pts_i[100][2], center_i[100][2];
    int nptsi=1;
    int resol_i[100];
    string connectors_i[100];
    
    double dwall 	  = (*Machine).inboard_generic_build(2);  // IB blanket thickness
	double dhtrepvert = (*Machine).vertical_generic_build(3); // repl HT shield thickness
	double newtopref  = etahiplsm + 1.2577; // to replace OB blanket ref (etahi2)


// ... Define the Outer Contour of Top Shield ...
// ... ... reference points ... ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 9;
     double xs[10], ys[10];
	 
	// 1 Bottom-left corner (outside, goes CW)
      xs[0] = ksilo1; 
      ys[0] = etahi1;

	// 2 Top-left, start of slope      
      xs[1] = xs[0]; 
      ys[1] = newtopref - dhtrepvert; // etahi2 - dhtrepvert; // replaced - 0.5*dwall
      // need to ref top of plasma + 1.2577 factor taken from OB FW part description
      // since repl HT can no longer ref OB blanket due to manifolding, mod 1/30/12
      // top of FW = top of plasma + 1.2577 = bottom of repl HT
	
	// 3 Middle of slope, determines angle      
      xs[2] = xs[1] + (0.5/2.0) *dwall;
      ys[2] = newtopref - (0.5/2.0) *dwall; // etahi2 - (0.5/2.0) *dwall; 
	
	// 4 Upper top-left corner, ending slope      
      xs[3] = xs[1] + (0.5) *dwall; // removed + 1/3 
      ys[3] = newtopref; // etahi2;

/*cout << "xs[0], ys[0] = " << xs[0] << ",	" << ys[0] << "\n";
cout << "xs[1], ys[1] = " << xs[1] << ",	" << ys[1] << "\n";
cout << "xs[2], ys[2] = " << xs[2] << ",	" << ys[2] << "\n";
cout << "xs[3], ys[3] = " << xs[3] << ",	" << ys[3] << "\n";    
*/
	// 5 top-right corner, invisible extention to push out HT shield for manifolding
      xs[4] = ksilo2+0.7; // arbitrary for now
      ys[4] = ys[3]; // etahi2;

	// 6 added - top-right correct corner, invisible extention to push out HT shield 
      xs[5] = ksilo2;
      ys[5] = ys[3];

	// 7 bottom-right corner of upper shield
      xs[6] = ksilo2;
      ys[6] = newtopref - dhtrepvert; // etahi2 - dhtrepvert; // removed - 0.5*dwall

	// 8 top-left inside 90 deg corner      
      xs[7] = ksilo1 + dhtrepvert; // removed - 0.333*dwall
      ys[7] = newtopref - dhtrepvert; // etahi2 - dhtrepvert; // removed - 0.5*dwall

	// 9 end on top of IB blanket      
      xs[8] = xs[7];
      ys[8] = ys[0];


     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
      if(i == 1 || i == 5) connectors_o[i] = "line"; } // || means or (is this line redundant)

// ... Update the Top Shield ...
	(*Shield)[0].get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, 
										  ref_pts_o, connectors_i, connectors_o);

// Define the Inner Target Curve For Blanket Shape
	double targcurve[10000][2];  int tgl = 0;

// ... Connect the Reference Points into Contours ...
	(*Shield)[0].connect(center_i, center_o, targcurve, tgl);


// ... Establish Contour Depth for volume calculation ...
	const double pi = 4.0*atan(1.0); // pi = 3.14159....
	double * d_in; double * d_out;
	int lin  = (*Shield)[0].contour_size("inner"); // 0
	int lout = (*Shield)[0].contour_size("outer"); // 160-180

	d_in  = new double [lin];
	d_out = new double [lout];
	
 	for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 	for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Shield)[0].contour(i, 0, "outer")/16.0;
 	
// ... Volume and Surface Areas for Top Shield
	(*Shield)[0].get_depth(d_in, d_out);
	(*Shield)[0].get_volume();
	(*Shield)[0].get_surface_areas();


//--------------------
// Bottom Shield
// ... Define the Contours of Bottom Shield ...
(*Shield)[1].get_contour_size(lin,  "inner");
(*Shield)[1].get_contour_size(lout, "outer");
 double (*ctr_in)[2], (*ctr_out)[2];
 ctr_in  = new double [lin][2];
 ctr_out = new double [lout][2];

for(int i=0; i<lin; i++)
{ ctr_in[i][0] = 0.0;
  ctr_in[i][1] = 0.0;
  d_in[i] 	   = 0.0; }

for(int i=0; i<lout; i++)
{  ctr_out[i][0] =  (*Shield)[0].contour(lout-1-i, 0, "outer");
   ctr_out[i][1] = -(*Shield)[0].contour(lout-1-i, 1, "outer");
   d_out[i] 	 = 2.0*pi*ctr_out[i][0]/16.0; }
(*Shield)[1].get_contours(ctr_in, ctr_out, d_in, d_out);


// ... Volume and Surface Areas for bottom shield ............
(*Shield)[1].get_volume();
(*Shield)[1].get_surface_areas();

delete [] d_in;
delete [] d_out;
delete [] ctr_in;
delete [] ctr_out;

} // eof


// ====================================================================================
// HT Shield - SCLL - splits into IB and OB but builds it continuously

void design_ht_shield_sic(DesignPoint * Machine, 	Part (*Shield)[2], 	Part * RefLeft,
                      			Part (*RefDiag)[2], Part * RefTop, 		Part * RefRight)
                      			// RefLeft is		BlanketInboard
                      			// RefDiag[2] is	RepHTShield, [0] top, [1] bottom
                      			// RefTop is		BlanketOutboard
                      			// RefRight is 		BlanketII
{ 
if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour of HT Shield...
    double ksilo1 = (*RefLeft).bounding_box(0, "outer");    // LHS of IB blanket
    double ksilo2 = (*RefRight).bounding_box(0, "outer");   // IB of blanket II near top
    double ksihi  = (*RefRight).bounding_box(2, "outer");   // OB of blanket II
    double etahi  = (*RefRight).bounding_box(3, "outer");   // top of blanket II (y)
    double etahi2 = (*RefDiag)[0].bounding_box(3, "outer"); // top of repl HT shield (y)
    double ksihi2 = (*RefDiag)[0].bounding_box(2, "outer"); // RHS of repl HT shield (x)
//cout << "ksilo1 = " << ksilo1 << "\n";
//cout << "ksihi  = " << ksihi << "\n";
//cout << (*RefDiag)[0].bounding_box(0, "outer") << "\n";


	double gapinb 	  = (*Machine).inboard_generic_build(7);   // IB blanket and HT shield gap
    double dhtinb 	  = (*Machine).inboard_generic_build(8);   // IB HT shield thickness, modified to be a little less in Aries.cpp	
	double dbwoutb2	  = (*Machine).outboard_generic_build(2);  // OB blanket I  thickness
	double dbwoutb 	  = (*Machine).outboard_generic_build(8);  // OB blanket II thickness
	double gapoutb    = (*Machine).outboard_generic_build(12); // blanket II to HT shield gap
    double dhtoutb    = (*Machine).outboard_generic_build(13); // OB HT shield thickness, modified little less in Aries.cpp
	double dhtrepvert = (*Machine).vertical_generic_build(3);  // repl HT shield thickness
	double gapvert 	  = (*Machine).vertical_generic_build(4);  // repl HT shield and HT shield gap
	double dhtvert 	  = (*Machine).vertical_generic_build(5);  // vertical HT shield thickness

// ... reference points ...
     double xs[4], ys[4];
     double ref_pts_i[100][2];
     int nptsi = 4;

	// 1 Start inner contour at IB midplane
     xs[0] = ksilo1 - gapinb; // LHS of IB blanket - gap
     ys[0] = 0.0;

	// 2 Top-left start of slope     
     xs[1] = xs[0]; // removed ksilo2 + 1.875*dhtinb + gapinb; 
     ys[1] = etahi2 - dhtrepvert + gapinb; // top of repl HT - repl HT thickness + gap // etahi - dhtrepvert+0.4;

	// 3 Top-right of slope     
     xs[2] = ksilo1 + dhtrepvert; // 3.665; // want 5.7
     ys[2] = etahi2 + gapvert; // 3.3; // want 4.6
     
	// 2 Top-left corner on top of repl HT shield - mod 1/31/12 for manifolding area 
//     xs[1] = xs[0]; // removed ksilo2 + 1.875*dhtinb + gapinb; 
//     ys[1] = etahi2 + gapvert; // removed - dhtrepvert; 
     
	// 3 Top-right of inner HT shield
	// LCC added top-right of inner HT shield
//   xs[2] = ksilo2 + dbwoutb + gapinb; 
//   ys[2] = etahi + gapvert;

	// 4 Bottom-right end at mid-plane     
     xs[3] = ksihi + gapinb; 
     ys[3] = 0.0;

/*cout << "HT Shield Inner: \n";
cout << "xs[0], ys[0] = " << xs[0] << ",	" << ys[0] << "\n";
cout << "xs[1], ys[1] = " << xs[1] << ",	" << ys[1] << "\n";
cout << "xs[2], ys[2] = " << xs[2] << ",	" << ys[2] << "\n";
cout << "xs[3], ys[3] = " << xs[3] << ",	" << ys[3] << "\n\n";
*/

     for(int i=0; i<nptsi; i++)
     { ref_pts_i[i][0] = xs[i];
       ref_pts_i[i][1] = ys[i]; }

// Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptsi; i++)
      	{ ref_pts_i[nptsi+i][0] =  ref_pts_i[nptsi-i-2][0];
          ref_pts_i[nptsi+i][1] = -ref_pts_i[nptsi-i-2][1]; }
          nptsi = 2*nptsi-1;

// Patch up the last segment with a line
       ref_pts_i[nptsi-1][0] = ref_pts_i[0][0];
       ref_pts_i[nptsi-1][1] = ref_pts_i[0][1] - 0.1; // was 0.5

// ... resolution ...
    int resol_i[100];
    for(int i=0; i< nptsi; i++)
    { resol_i[i] = 20; } // 20

// ... connectors ...
    string connectors_i[100];
    for(int i=0; i<nptsi; i++)
    { connectors_i[i]       = "target"; }
      connectors_i[nptsi-1] = "line"; // last patched segment is a line


// ----------------------------
// ... Outer Loop ...
// ... reference points ...
     double ref_pts_o[100][2];
     int nptso=6;

/* dcll here for reference
     ref_pts_o[1][0] = ref_pts_i[0][0] - dhtinb;
     ref_pts_o[1][1] = etahi - dhtrepvert + 0.5*dhtvert;
     
     ref_pts_o[2][0] = ksilo1 + dhtrepvert;
     ref_pts_o[2][1] = ref_pts_i[2][1] + dhtvert;
     
     ref_pts_o[3][0] = ref_pts_i[2][0];
     ref_pts_o[3][1] = ref_pts_i[2][1] + dhtvert;
     
     ref_pts_o[4][0] = ref_pts_i[2][0] + dhtoutb;
     ref_pts_o[4][1] = ref_pts_i[2][1];
     
     ref_pts_o[5][0] = ref_pts_i[3][0] + dhtoutb;
     ref_pts_o[5][1] = 0.0;
*/

	// 1 Bottom-left at midplane, ok
	 ref_pts_o[0][0] = xs[0] - dhtinb; // xs[0] - IB HT shield thickness
     ref_pts_o[0][1] = 0.0;

	// 2 Top-left, start of repl HT shield slope 
     ref_pts_o[1][0] = xs[0] - dhtinb; 
     ref_pts_o[1][1] = etahi2; // + dhtvert/2.0; // - dhtinb; // top of repl HT (y) + thickness

	// 3 Top-left, end of slope, mod 1/31/12 for manifolding area
	// would like to remove this and make the LH corner sharp but get a bus error...
     //ref_pts_o[2][0] = ksilo1 + dhtrepvert; // removed + dhtinb; 
     //ref_pts_o[2][1] = ref_pts_i[1][1] + dhtrepvert + gapvert + dhtvert; // removed + 2.0*dhtrepvert;
	 ref_pts_o[2][0] = xs[0] + dhtinb + gapinb; // was ksilo1 + dhtrepvert; 
     ref_pts_o[2][1] = etahi2 + dhtvert + gapvert; // + dhtrepvert + gapvert
	
	// 4 Top-right of HT shield
     ref_pts_o[3][0] = ksihi2 + gapoutb + dhtoutb+0.1; // extend out for manifolding, ksilo2 + dbwoutb + gapoutb + 0.1; // was ksilo2 + dbwoutb + gapoutb;  // replaced ref_pts_i[1][0]; // not 2,0
     ref_pts_o[3][1] = etahi2 + dhtvert + gapvert; // arbitrary ref_pts_o[2][1]+0.2; // ref_pts_i[1][1] + dhtvert + 2.0*dhtrepvert; // not 2,1

	// 5 Blanket II inside corner ref, removed for mod 1/31/12 for manifolding area 
     ref_pts_o[4][0] = ref_pts_o[3][0]; // ksilo2 + dbwoutb + gapoutb + dhtoutb;
     ref_pts_o[4][1] = etahi + dhtoutb/2.0; // ref_pts_i[1][1]; // top blanket II

	// 6 End at midplane
     ref_pts_o[5][0] = ksihi + gapoutb + dhtoutb; // ref_pts_i[2][0] + dhtoutb; // + dhtrepvert;
     ref_pts_o[5][1] = 0.0;

/*cout << "HT Shield Outer: \n";
cout << "ref_pts_o[0][0], ref_pts_o[0][1] = " << ref_pts_o[0][0] << ",	" << ref_pts_o[0][1] << "\n";
cout << "ref_pts_o[1][0], ref_pts_o[1][1] = " << ref_pts_o[1][0] << ",	" << ref_pts_o[1][1] << "\n";
cout << "ref_pts_o[2][0], ref_pts_o[2][1] = " << ref_pts_o[2][0] << ",	" << ref_pts_o[2][1] << "\n";
cout << "ref_pts_o[3][0], ref_pts_o[3][1] = " << ref_pts_o[3][0] << ",	" << ref_pts_o[3][1] << "\n";
cout << "ref_pts_o[4][0], ref_pts_o[4][1] = " << ref_pts_o[4][0] << ",	" << ref_pts_o[4][1] << "\n";
cout << "ref_pts_o[5][0], ref_pts_o[5][1] = " << ref_pts_o[5][0] << ",	" << ref_pts_o[5][1] << "\n";
*/


// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-1;
       

// Patch the last segment with a line
     ref_pts_o[nptso-1][0] = ref_pts_o[0][0];
     ref_pts_o[nptso-1][1] = ref_pts_o[0][1] - 0.1; // was 0.5

// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; } // # points between ref points

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i]       = "target"; }
      connectors_o[0]       = "line";
      connectors_o[nptso-1] = "line";


// ... Update the HT Shield ...
for(int i=0; i<2; i++) (*Shield)[i].get_contour_descriptions(nptsi, nptso, resol_i, resol_o,
                                 ref_pts_i, ref_pts_o, connectors_i, connectors_o);


// -----------------------------------------------
// Define the Inner Target Curve For Blanket Shape
// (HT shield attaches itself to the blanket and follows contour)
	// RefLeft is		BlanketInboard
    // RefDiag[2] is	RepHTShield, [0] top, [1] bottom
    // RefTop is		BlanketOutboard
    // RefRight is 		BlanketII
 
   double targcurve[10000][2];  int tgl = 0;

// 1) Get the outmost portion of RefLeft (BlanketInboard) contour ...
   for(int i=0; i<(*RefLeft).contour_size("outer"); i++)
   {
// ... outer normal to the curve ...
    double tagm = sqrt(pow(((*RefLeft).contour(i+1, 0, "outer")-
                            (*RefLeft).contour(i-1, 0, "outer")), 2) +
                       pow(((*RefLeft).contour(i+1, 1, "outer")-
                            (*RefLeft).contour(i-1, 1, "outer")), 2));
	double nx = - ((*RefLeft).contour(i+1, 1, "outer")-
                   (*RefLeft).contour(i-1, 1, "outer"))/tagm;
	double ny =   ((*RefLeft).contour(i+1, 0, "outer")-
                   (*RefLeft).contour(i-1, 0, "outer"))/tagm;
	
	if(i==0) { nx = 0.0; ny = 0.0; }
	if( nx == -1.0 || i == 0)
	{      targcurve[tgl][0] = (*RefLeft).contour(i, 0, "outer");
           targcurve[tgl][1] = (*RefLeft).contour(i, 1, "outer");
           tgl++; }
    if( ny == 1.0 ) break; }


// 2) Get the outmost portion of RefDiag[0] (RepHTShield top) contour ...
   for(int i=0; i<(*RefDiag)[0].contour_size("outer"); i++)
   {
// ... ... outer normal to the curve ... ...
    double tagm = sqrt(pow(((*RefDiag)[0].contour(i+1, 0, "outer")-
                            (*RefDiag)[0].contour(i-1, 0, "outer")), 2) +
                       pow(((*RefDiag)[0].contour(i+1, 1, "outer")-
                            (*RefDiag)[0].contour(i-1, 1, "outer")), 2));

	double nx = - ( (*RefDiag)[0].contour(i+1, 1, "outer")-
                    (*RefDiag)[0].contour(i-1, 1, "outer"))/tagm;
	double ny =   ( (*RefDiag)[0].contour(i+1, 0, "outer")-
                    (*RefDiag)[0].contour(i-1, 0, "outer"))/tagm;

	if(i==0) { nx = 0.0; ny = 0.0; }
	if( nx < 0.0 || ny == 1.0 || i==0)
	{      targcurve[tgl][0] = (*RefDiag)[0].contour(i, 0, "outer");
           targcurve[tgl][1] = (*RefDiag)[0].contour(i, 1, "outer");
           tgl++; }
	if(nx == 1) break; }


// 3) Get the outmost portion of RefTop (BlanketOutboard) contour ...
   for(int i=0; i<(*RefTop).contour_size("outer"); i++)
   {
// ... ... outer normal to the curve ... ...
    double tagm = sqrt( pow( ( (*RefTop).contour(i+1, 0, "outer")-
                               (*RefTop).contour(i-1, 0, "outer") ), 2) +
                        pow( ( (*RefTop).contour(i+1, 1, "outer")-
                               (*RefTop).contour(i-1, 1, "outer") ), 2) );

	double nx = - ( (*RefTop).contour(i+1, 1, "outer")-
                    (*RefTop).contour(i-1, 1, "outer") )/tagm;
	double ny =   ( (*RefTop).contour(i+1, 0, "outer")-
                    (*RefTop).contour(i-1, 0, "outer") )/tagm;

	if(i==0) { nx = 0.0; ny = 0.0; }
	if( ny == 1.0)
	{      targcurve[tgl][0] = (*RefTop).contour(i, 0, "outer");
           targcurve[tgl][1] = (*RefTop).contour(i, 1, "outer");
           tgl++; } }

//cout << (*RefRight).contour_size("outer") << "\n"; // 493


// 4) Get the outmost portion of RefRight (BlanketII) contour ...
   for(int i=0; i<(*RefRight).contour_size("outer"); i++)
   {
// ... ... outer normal to the curve ... ...
    double tagm = sqrt( pow( ( (*RefRight).contour(i+1, 0, "outer")-
                               (*RefRight).contour(i-1, 0, "outer") ), 2) +
                        pow( ( (*RefRight).contour(i+1, 1, "outer")-
                               (*RefRight).contour(i-1, 1, "outer") ), 2) );

	double nx = - ( (*RefRight).contour(i+1, 1, "outer")-
                    (*RefRight).contour(i-1, 1, "outer") )/tagm;
	double ny =   ( (*RefRight).contour(i+1, 0, "outer")-
                    (*RefRight).contour(i-1, 0, "outer") )/tagm;

	if(i==0) { nx = 0.0; ny = 0.0; }
	if( nx >= 0.0 && (*RefRight).contour(i, 1, "outer") >= 0.0 && i != 0)
	{      targcurve[tgl][0] = (*RefRight).contour(i, 0, "outer"); // offsets HT shield from blanket II
           targcurve[tgl][1] = (*RefRight).contour(i, 1, "outer");
           tgl++; } }


// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<tgl; i++)
      { targcurve[tgl+i][0] =  targcurve[tgl-i-2][0];
        targcurve[tgl+i][1] = -targcurve[tgl-i-2][1]; }
       tgl = 2*tgl-2;


// ... Connect the Reference Points into Contours ...
double center_i[100][2], center_o[100][2];
for(int i=0; i<2; i++) (*Shield)[i].connect(center_i, center_o, targcurve, tgl);


// ... Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin  = (*Shield)[0].contour_size("inner");
 int lout = (*Shield)[0].contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

// LCC shield[0] and [1] are the same...
 for(int k=0; k<2; k++)
 {  for(int i=0; i<lin; i++)
    {  d_in[i] = 2.0*pi*(*Shield)[k].contour(i, 0, "inner")/16.0;
       if( k == 0 &&    (*Shield)[k].contour(i, 0, "inner") >  ksilo1) d_in[i] = 0.0;
       if( k == 1 &&    (*Shield)[k].contour(i, 0, "inner") <= ksilo1) d_in[i] = 0.0; }
    for(int i=0; i<lout; i++)
    {  d_out[i] = 2.0*pi*(*Shield)[k].contour(i, 0, "outer")/16.0;
       if( k == 0 &&     (*Shield)[k].contour(i, 0, "outer") >  ksilo1) d_out[i] = 0.0;
       if( k == 1 &&     (*Shield)[k].contour(i, 0, "outer") <= ksilo1) d_out[i] = 0.0; }
    (*Shield)[k].get_depth(d_in, d_out); }


// ... Volume and Surface Areas ...
for(int i=0; i<2; i++) (*Shield)[i].get_volume();
for(int i=0; i<2; i++) (*Shield)[i].get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// HT Shield - DCLL

// Maintains one continuous HT Shield only in IB side, DCLL2 does not have HT Shield
void design_ht_shield_dcll(DesignPoint * Machine, 	Part * Shield, 		Part * RefLeft,
                      	   Part (*RefDiag)[2], 		Part * RefTop, 		Part * RefRight)
                      	   // RefLeft is 	BackWallInboard
                      	   // RefDiag[2] is RepHTShield, [0] top, [1] bottom
                      	   // RefTop is 	BlanketOutboard
                      	   // RefRight is 	BackWallOutboard, all seem OK when compared to SiC
{
if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour of HT Shield...
	double ksilo1 	  = (*RefLeft).bounding_box( 0, "outer");
	double ksilo2 	  = (*RefRight).bounding_box(0, "outer");
	double ksihi  	  = (*RefRight).bounding_box(2, "outer");
	double etahi  	  = (*RefRight).bounding_box(3, "outer");

	double gapinb 	  = (*Machine).inboard_generic_build(7);   // IB back wall and HT shield gap
	double dhtinb 	  = (*Machine).inboard_generic_build(8);   // IB HT shield thickness
	double dbwoutb 	  = (*Machine).outboard_generic_build(9);  // OB back wall II thickness
	double gapoutb 	  = (*Machine).outboard_generic_build(12); // OB back wall II and HT shield gap
	double dhtoutb 	  = (*Machine).outboard_generic_build(13); // OB HT shield thickness
	double dhtrepvert = (*Machine).vertical_generic_build(3);  // vert repl HT shield thickness
	double gapvert 	  = (*Machine).vertical_generic_build(4);  // vert repl HT shield and HT shield gap
	double dhtvert 	  = (*Machine).vertical_generic_build(5);  // vert HT shield thickness

    
// ... reference points ...
     double xs[4], ys[4];
     double ref_pts_i[100][2];
	 int nptsi=4;

     xs[0] = ksilo1 - gapinb;
     ys[0] = 0.0;

     xs[1] = xs[0];
     ys[1] = etahi - dhtrepvert;
     
     xs[2] = ksilo2 + dbwoutb + gapoutb;
     ys[2] = etahi  + gapvert;
     
     xs[3] = ksihi + gapoutb;
     ys[3] = 0.0;

     for(int i=0; i<nptsi; i++)
     { ref_pts_i[i][0] = xs[i];
       ref_pts_i[i][1] = ys[i]; }
/*cout << "xs[0], ys[0] = " << xs[0] << ",	" << ys[0] << "\n";
cout << "xs[1], ys[1] = " << xs[1] << ",	" << ys[1] << "\n";
cout << "xs[2], ys[2] = " << xs[2] << ",	" << ys[2] << "\n";
cout << "xs[3], ys[3] = " << xs[3] << ",	" << ys[3] << "\n";
*/

// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptsi; i++)
      { ref_pts_i[nptsi+i][0] =  ref_pts_i[nptsi-i-2][0];
        ref_pts_i[nptsi+i][1] = -ref_pts_i[nptsi-i-2][1]; }
       nptsi = 2*nptsi-1;


// ... ... resolution ... ...
    int resol_i[100];
    for(int i=0; i< nptsi; i++)
    { resol_i[i] = 40; }


// ... ... connectors ... ...
    string connectors_i[100];
    for(int i=0; i<nptsi; i++)
    { connectors_i[i] = "target"; }
      connectors_i[0] = "line";
      connectors_i[5] = "line";


// ... Outer Loop ... reference points ...
     double ref_pts_o[100][2];
     int nptso=6;

     ref_pts_o[0][0] = ref_pts_i[0][0] - dhtinb;
     ref_pts_o[0][1] = 0.0;
     
     ref_pts_o[1][0] = ref_pts_i[0][0] - dhtinb;
     ref_pts_o[1][1] = etahi - dhtrepvert + 0.5 * dhtvert;
     
     ref_pts_o[2][0] = ksilo1 + dhtrepvert;
     ref_pts_o[2][1] = ref_pts_i[2][1] + dhtvert;
     
     ref_pts_o[3][0] = ref_pts_i[2][0];
     ref_pts_o[3][1] = ref_pts_i[2][1] + dhtvert;
     
     ref_pts_o[4][0] = ref_pts_i[2][0] + dhtoutb;
     ref_pts_o[4][1] = ref_pts_i[2][1];
     
     ref_pts_o[5][0] = ref_pts_i[3][0] + dhtoutb;
     ref_pts_o[5][1] = 0.0;


// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-1;
       
// Patch the last segment with a line
     ref_pts_o[nptso-1][0] = ref_pts_o[0][0];
     ref_pts_o[nptso-1][1] = ref_pts_o[0][1] - 0.1;
    
// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }


// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] 		= "target"; }
    connectors_o[0] 		= "line";	connectors_o[1] 		= "line";
    connectors_o[2] 		= "line";	connectors_o[3] 		= "line";
    connectors_o[6] 		= "line";	connectors_o[7] 		= "line";
    connectors_o[8] 		= "line";	connectors_o[9] 		= "line";
    connectors_o[nptso-1] 	= "line";


// ... Update the HT Shield ...
 (*Shield).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                    connectors_i, connectors_o);

// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2]; int tgl = 0;

// ... get the outmost portion of RefLeft contour ...
   for(int i=0; i<(*RefLeft).contour_size("outer"); i++)
   {

// ... outer normal to the curve ...
    double tagm = sqrt(pow(((*RefLeft).contour(i+1, 0, "outer") -
                            (*RefLeft).contour(i-1, 0, "outer")), 2) +
                       pow(((*RefLeft).contour(i+1, 1, "outer") -
                            (*RefLeft).contour(i-1, 1, "outer")), 2));
	double nx 	= 		 - ((*RefLeft).contour(i+1, 1, "outer")-
                   			(*RefLeft).contour(i-1, 1, "outer"))/tagm;
	double ny 	=   	   ((*RefLeft).contour(i+1, 0, "outer")-
                   			(*RefLeft).contour(i-1, 0, "outer"))/tagm;
	
	if(i==0) { nx = 0.0; ny = 0.0; }
	if( nx == -1.0 || i == 0)
	{      targcurve[tgl][0] = (*RefLeft).contour(i, 0, "outer");
           targcurve[tgl][1] = (*RefLeft).contour(i, 1, "outer");
           tgl++; }
      if( ny == 1.0 ) break; }


// ... get the outmost portion of RefDiag[0] contour ...
   for(int i=0; i<(*RefDiag)[0].contour_size("outer"); i++)
   {

// ... outer normal to the curve ...
    double tagm = sqrt( pow( ( (*RefDiag)[0].contour(i+1, 0, "outer")-
                               (*RefDiag)[0].contour(i-1, 0, "outer") ), 2) +
                        pow( ( (*RefDiag)[0].contour(i+1, 1, "outer")-
                               (*RefDiag)[0].contour(i-1, 1, "outer") ), 2) );
	double nx 	=  		   - ( (*RefDiag)[0].contour(i+1, 1, "outer")-
                   			   (*RefDiag)[0].contour(i-1, 1, "outer") )/tagm;
	double ny 	=   	     ( (*RefDiag)[0].contour(i+1, 0, "outer")-
                      		   (*RefDiag)[0].contour(i-1, 0, "outer") )/tagm;

	if(i==0) { nx = 0.0; ny = 0.0; }
	if( nx < 0.0 || ny == 1.0 || i==0)
	{      targcurve[tgl][0] = (*RefDiag)[0].contour(i, 0, "outer");
           targcurve[tgl][1] = (*RefDiag)[0].contour(i, 1, "outer");
           tgl++; }
	if(nx == 1) break; }


// ... get the outmost portion of RefTop contour ...
   for(int i=0; i<(*RefTop).contour_size("outer"); i++)
   {

// ... outer normal to the curve ...
    double tagm = sqrt( pow( ( (*RefTop).contour(i+1, 0, "outer")-
                               (*RefTop).contour(i-1, 0, "outer") ), 2) +
                        pow( ( (*RefTop).contour(i+1, 1, "outer")-
                               (*RefTop).contour(i-1, 1, "outer") ), 2) );
	double nx 	=  		   - ( (*RefTop).contour(i+1, 1, "outer")-
    	                       (*RefTop).contour(i-1, 1, "outer") )/tagm;
	double ny 	= 		     ( (*RefTop).contour(i+1, 0, "outer")-
                      		   (*RefTop).contour(i-1, 0, "outer") )/tagm;

	if(i==0) { nx = 0.0; ny = 0.0; }
	if( ny == 1.0)
	{      targcurve[tgl][0] = (*RefTop).contour(i, 0, "outer");
           targcurve[tgl][1] = (*RefTop).contour(i, 1, "outer");
           tgl++; } }


// ... get the outmost portion of RefRight contour ...
   for(int i=0; i<(*RefRight).contour_size("outer"); i++)
   {

// ... outer normal to the curve ...
    double tagm = sqrt( pow( ( (*RefRight).contour(i+1, 0, "outer")-
                               (*RefRight).contour(i-1, 0, "outer") ), 2) +
                        pow( ( (*RefRight).contour(i+1, 1, "outer")-
                               (*RefRight).contour(i-1, 1, "outer") ), 2) );
	double nx 	= 		   - ( (*RefRight).contour(i+1, 1, "outer")-
                    		   (*RefRight).contour(i-1, 1, "outer") )/tagm;
	double ny 	= 	 		 ( (*RefRight).contour(i+1, 0, "outer")-
                    		   (*RefRight).contour(i-1, 0, "outer") )/tagm;

	if(i==0) { nx = 0.0; ny = 0.0; }
	if( nx >= 0.0 && (*RefRight).contour(i, 1, "outer") >= 0.0 && i != 0)
	{      targcurve[tgl][0] = (*RefRight).contour(i, 0, "outer");
           targcurve[tgl][1] = (*RefRight).contour(i, 1, "outer");
           tgl++; } }


// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<tgl; i++)
      { targcurve[tgl+i][0] =  targcurve[tgl-i-2][0];
        targcurve[tgl+i][1] = -targcurve[tgl-i-2][1]; }
       tgl = 2*tgl-2;

// ... Connect the Reference Points into Contours ...
 double center_i[100][2], center_o[100][2];
 (*Shield).connect(center_i, center_o, targcurve, tgl);


// ... Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin  = (*Shield).contour_size("inner");
 int lout = (*Shield).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

    for(int i=0; i<lin; i++)
    	{ d_in[i]  = 2.0*pi*(*Shield).contour(i, 0, "inner")/16.0; 
    	//cout << d_in[i] << "\n"; 
    	}
    for(int i=0; i<lout; i++)
    	{ d_out[i] = 2.0*pi*(*Shield).contour(i, 0, "outer")/16.0; 
    	//cout << d_out[i] << "\n"; 
    	}
    	
    (*Shield).get_depth(d_in, d_out);

// ... Volume and Surface Areas ...
 (*Shield).get_volume();
 (*Shield).get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Skeleton Ring - DCLL2

// For DCLL2, make the RefBody of SkeletonRing that of HT Shield DCLL
// Essentially turn HT Shield into SkeletonRing (No HT Shield for DCLL2)

void design_skeleton_ring2(DesignPoint * Machine, 	Part * Ring, 		Part * RefLeft,
                      	   Part (*RefDiag)[2], 		Part * RefTop, 		Part * RefRight)
                      	   // RefLeft is 	BackWallInboard
                      	   // RefDiag[2] is RepHTShield, [0] top, [1] bottom
                      	   // RefTop is 	BlanketOutboard
                      	   // RefRight is 	BackWallOutboard
{
if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour of HT Shield...
	double ksilo1 	  = (*RefLeft).bounding_box( 0, "outer"); // BW IB edge
	double ksilo2 	  = (*RefRight).bounding_box(0, "outer"); // BW OB ?
	double ksihi  	  = (*RefRight).bounding_box(2, "outer"); // BW OB ?
	double etahi  	  = (*RefRight).bounding_box(3, "outer"); // BW OB top

// from HT Shield function
	double gapinb 	  = (*Machine).inboard_generic_build(9);   // IB back wall and SkelRing gap
	double dsrinb 	  = (*Machine).inboard_generic_build(10);  // IB SkelRing thickness
	double dbwoutb 	  = (*Machine).outboard_generic_build(9);  // OB back wall II thickness
	double gapoutb 	  = (*Machine).outboard_generic_build(14); // OB back wall II and SkelRing gap
	double dsroutb 	  = (*Machine).outboard_generic_build(15); // OB SkelRing thickness
	double drepvert   = (*Machine).vertical_generic_build(3);  // vert repl HT shield thickness
	double gapvert 	  = (*Machine).vertical_generic_build(6);  // vert repl HT shield and SkelRing gap
	double dsrvert 	  = (*Machine).vertical_generic_build(7);  // vert SkelRing thickness
  
// ... reference points ...
     double xs[4], ys[4];
     double ref_pts_i[100][2];
	 int nptsi=4;

     xs[0] = ksilo1 - gapinb;
     ys[0] = 0.0;

     xs[1] = xs[0];
     ys[1] = etahi - drepvert;
     
     xs[2] = ksilo2 + dbwoutb + gapoutb;
     ys[2] = etahi  + gapvert;
     
     xs[3] = ksihi + gapoutb;
     ys[3] = 0.0;

     for(int i=0; i<nptsi; i++)
     { ref_pts_i[i][0] = xs[i];
       ref_pts_i[i][1] = ys[i]; }
       
/*cout << "xs[0], ys[0] = " << xs[0] << ",	" << ys[0] << "\n";
cout << "xs[1], ys[1] = " << xs[1] << ",	" << ys[1] << "\n";
cout << "xs[2], ys[2] = " << xs[2] << ",	" << ys[2] << "\n";
cout << "xs[3], ys[3] = " << xs[3] << ",	" << ys[3] << "\n";
*/

// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptsi; i++)
      { ref_pts_i[nptsi+i][0] =  ref_pts_i[nptsi-i-2][0];
        ref_pts_i[nptsi+i][1] = -ref_pts_i[nptsi-i-2][1]; }
       nptsi = 2*nptsi-1;


// ... ... resolution ... ...
    int resol_i[100];
    for(int i=0; i< nptsi; i++)
    { resol_i[i] = 40; }


// ... ... connectors ... ...
    string connectors_i[100];
    for(int i=0; i<nptsi; i++)
    { connectors_i[i] = "target"; }
      connectors_i[0] = "line";
      connectors_i[5] = "line";


// ... Outer Loop ... reference points ...
     double ref_pts_o[100][2];
     int nptso=6;

     ref_pts_o[0][0] = ref_pts_i[0][0] - dsrinb;
     ref_pts_o[0][1] = 0.0;
     
     ref_pts_o[1][0] = ref_pts_i[0][0] - dsrinb;
     ref_pts_o[1][1] = etahi - drepvert + 0.5*dsrvert;
     
     ref_pts_o[2][0] = ksilo1 + drepvert;
     ref_pts_o[2][1] = ref_pts_i[2][1] + dsrvert;
     
     ref_pts_o[3][0] = ref_pts_i[2][0];
     ref_pts_o[3][1] = ref_pts_i[2][1] + dsrvert;
     
     ref_pts_o[4][0] = ref_pts_i[2][0] + dsroutb;
     ref_pts_o[4][1] = ref_pts_i[2][1];
     
     ref_pts_o[5][0] = ref_pts_i[3][0] + dsroutb;
     ref_pts_o[5][1] = 0.0;


// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-1;
       
// Patch the last segment with a line
     ref_pts_o[nptso-1][0] = ref_pts_o[0][0];
     ref_pts_o[nptso-1][1] = ref_pts_o[0][1] - 0.1;
    
// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }


// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] 		= "target"; }
    connectors_o[0] 		= "line";	connectors_o[1] 		= "line";
    connectors_o[2] 		= "line";	connectors_o[3] 		= "line";
    connectors_o[6] 		= "line";	connectors_o[7] 		= "line";
    connectors_o[8] 		= "line";	connectors_o[9] 		= "line";
    connectors_o[nptso-1] 	= "line";


// ... Update the Skeleton Ring ...
 (*Ring).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                    connectors_i, connectors_o);

// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2]; int tgl = 0;

// ... get the outmost portion of RefLeft contour ...
   for(int i=0; i<(*RefLeft).contour_size("outer"); i++)
   {

// ... outer normal to the curve ...
    double tagm = sqrt(pow(((*RefLeft).contour(i+1, 0, "outer") -
                            (*RefLeft).contour(i-1, 0, "outer")), 2) +
                       pow(((*RefLeft).contour(i+1, 1, "outer") -
                            (*RefLeft).contour(i-1, 1, "outer")), 2));
	double nx 	= 		 - ((*RefLeft).contour(i+1, 1, "outer")-
                   			(*RefLeft).contour(i-1, 1, "outer"))/tagm;
	double ny 	=   	   ((*RefLeft).contour(i+1, 0, "outer")-
                   			(*RefLeft).contour(i-1, 0, "outer"))/tagm;
	
	if(i==0) { nx = 0.0; ny = 0.0; }
	if( nx == -1.0 || i == 0)
	{      targcurve[tgl][0] = (*RefLeft).contour(i, 0, "outer");
           targcurve[tgl][1] = (*RefLeft).contour(i, 1, "outer");
           tgl++; }
      if( ny == 1.0 ) break; }


// ... get the outmost portion of RefDiag[0] contour ...
   for(int i=0; i<(*RefDiag)[0].contour_size("outer"); i++)
   {

// ... outer normal to the curve ...
    double tagm = sqrt( pow( ( (*RefDiag)[0].contour(i+1, 0, "outer")-
                               (*RefDiag)[0].contour(i-1, 0, "outer") ), 2) +
                        pow( ( (*RefDiag)[0].contour(i+1, 1, "outer")-
                               (*RefDiag)[0].contour(i-1, 1, "outer") ), 2) );
	double nx 	=  		   - ( (*RefDiag)[0].contour(i+1, 1, "outer")-
                   			   (*RefDiag)[0].contour(i-1, 1, "outer") )/tagm;
	double ny 	=   	     ( (*RefDiag)[0].contour(i+1, 0, "outer")-
                      		   (*RefDiag)[0].contour(i-1, 0, "outer") )/tagm;

	if(i==0) { nx = 0.0; ny = 0.0; }
	if( nx < 0.0 || ny == 1.0 || i==0)
	{      targcurve[tgl][0] = (*RefDiag)[0].contour(i, 0, "outer");
           targcurve[tgl][1] = (*RefDiag)[0].contour(i, 1, "outer");
           tgl++; }
	if(nx == 1) break; }


// ... get the outmost portion of RefTop contour ...
   for(int i=0; i<(*RefTop).contour_size("outer"); i++)
   {

// ... outer normal to the curve ...
    double tagm = sqrt( pow( ( (*RefTop).contour(i+1, 0, "outer")-
                               (*RefTop).contour(i-1, 0, "outer") ), 2) +
                        pow( ( (*RefTop).contour(i+1, 1, "outer")-
                               (*RefTop).contour(i-1, 1, "outer") ), 2) );
	double nx 	=  		   - ( (*RefTop).contour(i+1, 1, "outer")-
    	                       (*RefTop).contour(i-1, 1, "outer") )/tagm;
	double ny 	= 		     ( (*RefTop).contour(i+1, 0, "outer")-
                      		   (*RefTop).contour(i-1, 0, "outer") )/tagm;

	if(i==0) { nx = 0.0; ny = 0.0; }
	if( ny == 1.0)
	{      targcurve[tgl][0] = (*RefTop).contour(i, 0, "outer");
           targcurve[tgl][1] = (*RefTop).contour(i, 1, "outer");
           tgl++; } }


// ... get the outmost portion of RefRight contour ...
   for(int i=0; i<(*RefRight).contour_size("outer"); i++)
   {

// ... outer normal to the curve ...
    double tagm = sqrt( pow( ( (*RefRight).contour(i+1, 0, "outer")-
                               (*RefRight).contour(i-1, 0, "outer") ), 2) +
                        pow( ( (*RefRight).contour(i+1, 1, "outer")-
                               (*RefRight).contour(i-1, 1, "outer") ), 2) );
	double nx 	= 		   - ( (*RefRight).contour(i+1, 1, "outer")-
                    		   (*RefRight).contour(i-1, 1, "outer") )/tagm;
	double ny 	= 	 		 ( (*RefRight).contour(i+1, 0, "outer")-
                    		   (*RefRight).contour(i-1, 0, "outer") )/tagm;

	if(i==0) { nx = 0.0; ny = 0.0; }
	if( nx >= 0.0 && (*RefRight).contour(i, 1, "outer") >= 0.0 && i != 0)
	{      targcurve[tgl][0] = (*RefRight).contour(i, 0, "outer");
           targcurve[tgl][1] = (*RefRight).contour(i, 1, "outer");
           tgl++; } }


// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<tgl; i++)
      { targcurve[tgl+i][0] =  targcurve[tgl-i-2][0];
        targcurve[tgl+i][1] = -targcurve[tgl-i-2][1]; }
       tgl = 2*tgl-2;

// ... Connect the Reference Points into Contours ...
 double center_i[100][2], center_o[100][2];
 (*Ring).connect(center_i, center_o, targcurve, tgl);


// ... Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin  = (*Ring).contour_size("inner");
 int lout = (*Ring).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

    for(int i=0; i<lin; i++)
    	{ d_in[i]  = 2.0*pi*(*Ring).contour(i, 0, "inner")/16.0; 
    	//cout << d_in[i] << "\n"; 
    	}
    for(int i=0; i<lout; i++)
    	{ d_out[i] = 2.0*pi*(*Ring).contour(i, 0, "outer")/16.0; 
    	//cout << d_out[i] << "\n"; 
    	}
    	
    (*Ring).get_depth(d_in, d_out);

// ... Volume and Surface Areas ...
 (*Ring).get_volume();
 (*Ring).get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Skeleton Ring - DCLL

// NOTE this assumes HT shield (HTShieldDCLL) is in one piece
void design_skeleton_ring(DesignPoint *Machine, Part *Ring, Part *RefBody) 
{ 												// RefBody is HTShieldDCLL
if( (*Machine).data_point_filter(0) > 0 ) return; 

// ... Define the Inner Contour of the Vacuum Vessel ...
// NEED to ensure these are correct, accounting for single HT shield body, not two
     double ksilo  = (*RefBody).bounding_box(0, "outer");
     double ksihi  = (*RefBody).bounding_box(2, "outer");
     double etahi  = (*RefBody).bounding_box(3, "outer");
     double ksimid = (ksilo + ksihi)/2.0;
     double reflen = (*RefBody).contour_size("outer");


// ... reference points ...
    double dhtinb     = (*Machine).inboard_generic_build(8);  // IB HT shield thickness  
	double gapinb     = (*Machine).inboard_generic_build(9);  // IB HT shield and SkelRing gap
	double dvinb      = (*Machine).inboard_generic_build(10); // IB SkelRing thickness
  //double gapoutb    = (*Machine).outboard_generic_build(14);// OB HT shield and SkelRing gap
    double dvoutb     = (*Machine).outboard_generic_build(15);// OB SkelRing thickness
    double drephtvert = (*Machine).vertical_generic_build(3); // vert repl HT shield thickness
    double dhtvert    = (*Machine).vertical_generic_build(5); // vert HT shield thickness
    double gapvert    = (*Machine).vertical_generic_build(6); // vert HT shield and SkelRing gap
	double dvvert     = (*Machine).vertical_generic_build(7); // vert SkelRing thickness


// ... reference points ...
      double ref_pts_i[100][2];
      double xs[5], ys[5];
      int nptsi=5;

      xs[0] = ksilo - gapinb;
      ys[0] = 0.0;

      xs[1] = xs[0];
      ys[1] = etahi - drephtvert - 0.5 * dhtvert; // 45 deg slope to mimic repl HT shield
      
      xs[2] = ksilo + dhtinb + drephtvert;
      ys[2] = etahi + gapvert;
      
      for(int i=0; i<reflen; i++)
      { if( (*RefBody).contour(i, 1, "outer") >= (etahi - dvvert) )
       { xs[3] = (*RefBody).contour(i, 0, "outer"); } }
         ys[3] = etahi + gapvert;
      
      xs[4] = ksihi;
      ys[4] = 0.0;
      
     for(int i=0; i<nptsi; i++)
     { ref_pts_i[i][0] = xs[i];
       ref_pts_i[i][1] = ys[i]; }


// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptsi; i++)
      { ref_pts_i[nptsi+i][0] =  ref_pts_i[nptsi-i-2][0];
        ref_pts_i[nptsi+i][1] = -ref_pts_i[nptsi-i-2][1]; }
      nptsi = 2*nptsi-1;

// ... ... resolution ... ...
    int resol_i[100];
    for(int i=0; i< nptsi; i++)
    { resol_i[i] = 20; }


// ... ... connectors ... ...
    string connectors_i[100];
    for(int i=0; i<nptsi; i++)
    { connectors_i[i] = "target"; }
    connectors_i[0] = "line";
    connectors_i[2] = "line";
    connectors_i[5] = "line";
    connectors_i[7] = "line";
    
    
// ... Outer Loop ... reference points ...
     double ref_pts_o[100][2];
     int nptso=6;

     ref_pts_o[0][0] = ref_pts_i[0][0] - dvinb;
     ref_pts_o[0][1] = 0.0;

     ref_pts_o[1][0] = ref_pts_o[0][0];
     ref_pts_o[1][1] = ref_pts_i[2][1] + 0.5 * dvvert;

     ref_pts_o[2][0] = ref_pts_i[2][0];
     ref_pts_o[2][1] = ref_pts_i[2][1] + dvvert;

     ref_pts_o[3][0] = ref_pts_i[3][0];
     ref_pts_o[3][1] = ref_pts_i[3][1] + dvvert;

     ref_pts_o[4][0] = ref_pts_i[3][0] + dvoutb;
     ref_pts_o[4][1] = ref_pts_i[3][1];
     
     ref_pts_o[5][0] = ref_pts_i[4][0] + dvoutb;
     ref_pts_o[5][1] = ref_pts_i[4][1];


// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-1;

// Patch the last segment with a line
     ref_pts_o[nptso-1][0] = ref_pts_o[0][0];
     ref_pts_o[nptso-1][1] = ref_pts_o[0][1] - 0.1;

// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
      if(i == 4 || i ==5) connectors_o[i] = "target"; }


// ... Update the Skeleton Ring ...
(*Ring).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, 
								 ref_pts_o, connectors_i, connectors_o);
                                 
                                 
// Define the Inner Target Curve For Skeleton Ring Shape
double targcurve[10000][2];
int tgl = (*RefBody).contour_size("outer");
  for(int i=0; i<tgl; i++)
  { targcurve[i][0] = (*RefBody).contour(i, 0, "outer");
    targcurve[i][1] = (*RefBody).contour(i, 1, "outer"); }


// Connect the Reference Points into Contours
double center_i[100][2], center_o[100][2];
(*Ring).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin  = (*Ring).contour_size("inner");
 int lout = (*Ring).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];


// Inboard loop of the skeleton ring
 for(int i=0; i<lin; i++)
 { d_in[i] = 2.0*pi*(*Ring).contour(i, 0, "inner")/16.0; }


// Outboard loop of the skeleton ring
 for(int i=0; i<lout; i++)
 { d_out[i] = 2.0*pi*(*Ring).contour(i, 0, "outer")/16.0; }
 (*Ring).get_depth(d_in, d_out);
 

// Volume and Surface Areas
for(int i=0; i<2; i++) (*Ring).get_volume();
for(int i=0; i<2; i++) (*Ring).get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Manifold - DCLL OB and vert in Laila's new (July 2011) design, contains PbLi+He
// Different than "manifolding" for routing of He & PbLi into/out of blankets and out of vessel

void design_manifold(DesignPoint * Machine, Part * Manifold, Part * RefBody) 
											// RefBody is skeleton ring
{
if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour of the Manifold ...
     double ksilo  = (*RefBody).bounding_box(0, "outer"); // IB of skel ring
     double ksihi  = (*RefBody).bounding_box(2, "outer"); // OB of skel ring
     double etahi  = (*RefBody).bounding_box(3, "outer"); // top (Y) of skel ring
     double reflen = (*RefBody).contour_size("outer");    // 466, length of contour
//cout << "ksilo = " << ksilo << "\n";
//cout << "ksihi = " << ksihi << "\n";
//cout << "etahi = " << etahi << "\n";

// ... ... reference points ... ...
     double ref_pts_i[100][2];
     int nptsi=5;
     double xs[5], ys[5];
     
     double dsrinb 	   = (*Machine).inboard_generic_build(10); // IB skel ring thickness
     double gapinb 	   = (*Machine).inboard_generic_build(11); // IB skel ring and manifold gap (SR was HT shield)
     double dvinb 	   = (*Machine).inboard_generic_build(12); // IB manifold thickness
     double gapoutb    = (*Machine).outboard_generic_build(16);// OB skel ring and manifold gap (SR was HT shield)
     double dvoutb 	   = (*Machine).outboard_generic_build(17);// OB manifold thickness
     double drephtvert = (*Machine).vertical_generic_build(3); // vert repl HT shield thickness
     double dsrvert    = (*Machine).vertical_generic_build(7); // vert skel ring thickness
     double gapvert    = (*Machine).vertical_generic_build(8); // vert skel ring and manifold gap (SR was HT shield)
     double dvvert 	   = (*Machine).vertical_generic_build(9); // vert manifold thickness


// ... reference points ...
      xs[0] = ksilo - gapinb; // start at IB midplane as if there is an IB manifold 
      ys[0] = 0.0;
      
      xs[1] = xs[0];  // bottom left of slope
      ys[1] = etahi - dsrvert - 0.025; // temp, removed - drephtvert - 0.5 * dhtvert
      
      xs[2] = ksilo + dsrinb + drephtvert; // top right of slope, removed + dhtinb + gapoutb
      ys[2] = etahi + gapvert;
      
      for(int i=0; i<reflen; i++)
    { if( (*RefBody).contour(i, 1, "outer") >= (etahi) ) // folloW SR contour, removed - dvvert
    { xs[3] = (*RefBody).contour(i, 0, "outer"); } } // top right of top, added - gapoutb
      ys[3] = etahi + gapvert;
      
      xs[4] = ksihi; // right midplane end
      ys[4] = 0.0;

/*cout << "xs[0], ys[0] = " << xs[0] << ", " << ys[0] << "\n";
cout << "xs[1], ys[1] = " << xs[1] << ", " << ys[1] << "\n";
cout << "xs[2], ys[2] = " << xs[2] << ", " << ys[2] << "\n";
cout << "xs[3], ys[3] = " << xs[3] << ", " << ys[3] << "\n";
cout << "xs[4], ys[4] = " << xs[4] << ", " << ys[4] << "\n";
*/
      
     for(int i=0; i<nptsi; i++)
     { ref_pts_i[i][0] = xs[i];
       ref_pts_i[i][1] = ys[i]; }


// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptsi; i++)
      { ref_pts_i[nptsi+i][0] =  ref_pts_i[nptsi-i-2][0];
        ref_pts_i[nptsi+i][1] = -ref_pts_i[nptsi-i-2][1]; }
        nptsi = 2*nptsi-1;


// ... ... resolution ... ...
    int resol_i[100];
    for(int i=0; i< nptsi; i++)
    	{ resol_i[i] = 20; }


// ... ... connectors ... ...
    string connectors_i[100];
    for(int i=0; i<nptsi; i++)
    { connectors_i[i] = "target"; }
      connectors_i[0] = "line";
      connectors_i[2] = "line";
      connectors_i[5] = "line";
      connectors_i[7] = "line";
    
    
// ... Outer Loop ...
// ... reference points ...
     double ref_pts_o[100][2];
     int nptso=6;

     ref_pts_o[0][0] = ref_pts_i[0][0] - dvinb;
     ref_pts_o[0][1] = 0.0;

//     ref_pts_o[1][0] = ref_pts_i[0][0] - dsrinb;
//     ref_pts_o[1][1] = etahi - drepvert + 0.5*dsrvert;

     ref_pts_o[1][0] = ref_pts_o[0][0];
     ref_pts_o[1][1] = ref_pts_i[2][1] - 0.5 * dsrvert; // removed + 0.5 * dvvert, helps slope stick to SR

     ref_pts_o[2][0] = ref_pts_i[2][0];
     ref_pts_o[2][1] = ref_pts_i[2][1] + dvvert;

     ref_pts_o[3][0] = ref_pts_i[3][0];
     ref_pts_o[3][1] = ref_pts_i[3][1] + dvvert;

     ref_pts_o[4][0] = ref_pts_i[3][0] + dvoutb;
     ref_pts_o[4][1] = ref_pts_i[3][1];
     
     ref_pts_o[5][0] = ref_pts_i[4][0] + dvoutb;
     ref_pts_o[5][1] = ref_pts_i[4][1];


/* ... Close the loop by adding a symmetric profile in the clockwise direction: ... */
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-1;

// Patch the last segment with a line
     ref_pts_o[nptso-1][0] = ref_pts_o[0][0];
     ref_pts_o[nptso-1][1] = ref_pts_o[0][1] - 0.1;


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }


// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
    if(i == 4 || i ==5) connectors_o[i] = "target"; }


// ... Update the Manifold ...
(*Manifold).get_contour_descriptions(nptsi, 	nptso, 		resol_i, 		resol_o,
									ref_pts_i,	ref_pts_o,	connectors_i,	connectors_o);
                 
                 
// Define the Inner Target Curve For Manifold Shape
double targcurve[10000][2];
int tgl = (*RefBody).contour_size("outer");
  for(int i=0; i<tgl; i++)
  { targcurve[i][0] = (*RefBody).contour(i, 0, "outer");
    targcurve[i][1] = (*RefBody).contour(i, 1, "outer"); }


// Connect the Reference Points into Contours
 double center_i[100][2], center_o[100][2];
(*Manifold).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0);
 double * d_in; double * d_out;
 int lin  = (*Manifold).contour_size("inner");;
 int lout = (*Manifold).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

// Inboard loop of the manifold
 for(int i=0; i<lin; i++)
 { d_in[i] = 2.0*pi*(*Manifold).contour(i, 0, "inner")/16.0; }

// Outboard loop of the manifold
 for(int i=0; i<lout; i++)
 { d_out[i] = 2.0*pi*(*Manifold).contour(i, 0, "outer")/16.0; }
 (*Manifold).get_depth(d_in, d_out);
 

// Volume and Surface Areas
 for(int i=0; i<2; i++) (*Manifold).get_volume();
 for(int i=0; i<2; i++) (*Manifold).get_surface_areas();
//cout << "(*Manifold).volume() = " << (*Manifold).volume() << "\n";

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Vacuum Vessel - SCLL  (single function)
// LCC - the "split" between "main" and "OB" VV occurs when it reduces from 40cm to 25cm

void design_vacuum_vessel(DesignPoint *Machine, Part (*Vessel)[4], Part (*RefBody)[2], Part (*Shield)[2])
 				// RefBody is HTShield; SiC uses Vessel[1-2] of [4], Shield is Repl HT shield
{ if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour of the Vacuum Vessel ...
     double ksilo  = (*RefBody)[0].bounding_box(0, "outer"); // LHS IB of HT shield
     double ksihi  = (*RefBody)[0].bounding_box(2, "outer"); // RHS OB of HT shield
     double etahi  = (*RefBody)[0].bounding_box(3, "outer"); // top (y) of HT shield
     double ksihi2 = (*Shield)[0].bounding_box(2, "outer");  // RHS of repl HT shield (x)
     double reflen = (*RefBody)[0].contour_size("outer");    // pts in contour (487)
     double ksimid = (ksilo + ksihi)/2.0; // middle of HT shield (?)
//cout << "ksilo  = " << ksilo << "\n";
//cout << "ksihi  = " << ksihi << "\n";
//cout << "etahi  = " << etahi << "\n";

// ... reference points ... from Aries.cpp
     double dwgap   = (*Machine).inboard_generic_build(13);  // IB HT shield and VV gap
     double dwall   = (*Machine).inboard_generic_build(14);  // IB VV thickness = 40cm
     double dwallv  = (*Machine).vertical_generic_build(11); // Vert VV thickness = 40cm
     double obhtthk = (*Machine).outboard_generic_build(13); // OB HT shield thickness
     double obvvgap = (*Machine).outboard_generic_build(18); // OB HT shield and VV gap
     double obvv    = (*Machine).outboard_generic_build(19); // OB VV thickness = 25cm

// ... reference points ...
      double xs[5], ys[5];
      double ref_pts_i[100][2];
      int nptsi = 4;

	// 1 Start at inner of IB VV at midplane   
      xs[0] = ksilo + dwgap; // Watch this (why + ?)
      ys[0] = 0.0;

	// 2 Top of VV near middle      
      xs[1] = (ksimid + ksilo)/2.0; // middle of HT shield + iner IB HT shield
      ys[1] = etahi + dwgap; // top of HT shield + gap

	// 3 Top-right of VV
	//  this if statement is making problems on the bottom...is it needed?
//	for(int i=0; i<reflen; i++)
//      { if( 	(*RefBody)[0].contour(i, 1, "outer") >= (etahi - dwall) ) // if HT shield is >= top of HT - IB VV thickness
//      { xs[2] = (*RefBody)[0].contour(i, 0, "outer") - 2.0*obhtthk; } } // OB HT thk // removed + dwall
//        ys[2] = etahi + dwgap;
	  xs[2] = ksihi2 + obvvgap + obhtthk + 0.125;
	  ys[2] = etahi + dwgap;	
	
	// 4 End at midplane       
      xs[3] = ksihi + obvvgap;
      ys[3] = 0.0;

/*cout << "VV inner: \n";
cout << "xs[0], ys[0] = " << xs[0] << ", " << ys[0] << "\n";
cout << "xs[1], ys[1] = " << xs[1] << ", " << ys[1] << "\n";
cout << "xs[2], ys[2] = " << xs[2] << ", " << ys[2] << "\n";
cout << "xs[3], ys[3] = " << xs[3] << ", " << ys[3] << "\n";
*/
      
     for(int i=0; i<nptsi; i++)
     { ref_pts_i[i][0] = xs[i];
       ref_pts_i[i][1] = ys[i]; }


// ... ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptsi; i++)
      { ref_pts_i[nptsi+i][0] =  ref_pts_i[nptsi-i-2][0];
        ref_pts_i[nptsi+i][1] = -ref_pts_i[nptsi-i-2][1]; }
        nptsi = 2*nptsi-1;


// Patch up the last segment with a line
       ref_pts_i[nptsi-1][0] = ref_pts_i[0][0];
       ref_pts_i[nptsi-1][1] = ref_pts_i[0][1] - 0.5;


// ... ... resolution ... ...
    int resol_i[100];
    for(int i=0; i< nptsi; i++)
    { resol_i[i] = 20; }


// ... ... connectors ... ...
    string connectors_i[100];
    for(int i=0; i<nptsi; i++)
    { connectors_i[i] 	  = "target"; }
    connectors_i[1] 	  = "line";
    connectors_i[4] 	  = "line";
    connectors_i[nptsi-1] = "line";
    
    
// ------------------------    
// ... Outer Loop ...
// ... ... reference points ... ...
     double ref_pts_o[100][2];
     int nptso=5;

	// 1 Start at midplane
     ref_pts_o[0][0] = ref_pts_i[0][0] - dwall;
     ref_pts_o[0][1] = 0.0;

	// 2 Top-left, use [1] ref point from inner contour
     ref_pts_o[1][0] = ref_pts_i[1][0];
     ref_pts_o[1][1] = ref_pts_i[1][1] + dwallv;

	// 3 Use [2] (top-right) ref point from inner contour
     ref_pts_o[2][0] = ref_pts_i[2][0];
     ref_pts_o[2][1] = ref_pts_i[2][1] + dwallv;
     
	// 4 Use [2] ref point from inner contour at
     ref_pts_o[3][0] = ref_pts_i[2][0] + obvv; // + (obvv+dwallv)/2.0; // thickness transition, removed + dwall*obvv/ibvv;
     ref_pts_o[3][1] = ref_pts_i[2][1];

	// 5 End at midplane
     ref_pts_o[4][0] = ref_pts_i[3][0] + obvv; // removed + dwall*obvv/ibvv;
     ref_pts_o[4][1] = 0.0;

/*cout << "VV outer: \n";
cout << "ref_pts_o[0][0], ref_pts_o[0][1] = " << ref_pts_o[0][0] << ", " << ref_pts_o[0][1] << "\n";
cout << "ref_pts_o[1][0], ref_pts_o[1][1] = " << ref_pts_o[1][0] << ", " << ref_pts_o[1][1] << "\n";
cout << "ref_pts_o[2][0], ref_pts_o[2][1] = " << ref_pts_o[2][0] << ", " << ref_pts_o[2][1] << "\n";
cout << "ref_pts_o[3][0], ref_pts_o[3][1] = " << ref_pts_o[3][0] << ", " << ref_pts_o[3][1] << "\n";
cout << "ref_pts_o[4][0], ref_pts_o[4][1] = " << ref_pts_o[4][0] << ", " << ref_pts_o[4][1] << "\n";
*/

// ... ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-1;

// Patch the last segment with a line
     ref_pts_o[nptso-1][0] = ref_pts_o[0][0];
     ref_pts_o[nptso-1][1] = ref_pts_o[0][1] - 0.5;

// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
      if(i == 0 || i == 3 || i ==4 || i == 7) connectors_o[i] = "target"; }


// ... Update the VV ...
for(int i=0; i<2; i++) (*Vessel)[i].get_contour_descriptions(nptsi, nptso, resol_i, 
		resol_o, ref_pts_i, ref_pts_o, connectors_i, connectors_o);
                                 
                                 
// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2];
int tgl = (*RefBody)[0].contour_size("outer");

  for(int i=0; i<tgl; i++)
  { targcurve[i][0] = (*RefBody)[0].contour(i, 0, "outer"); // RefBody is HT shield
    targcurve[i][1] = (*RefBody)[0].contour(i, 1, "outer"); }


// Connect the Reference Points into Contours
 double center_i[100][2], center_o[100][2];
 for(int i=0; i<2; i++) (*Vessel)[i].connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin  = (*Vessel)[0].contour_size("inner");;
 int lout = (*Vessel)[0].contour_size("outer");
 
 d_in  = new double [lin];
 d_out = new double [lout];


// Main body of the vacuum vessel
 for(int i=0; i<lin; i++)
 { d_in[i] = 2.0*pi*(*Vessel)[0].contour(i, 0, "inner")/16.0;
   if( (*Vessel)[0].contour(i, 0, "inner") > ref_pts_i[2][0] ) d_in[i] = 0.0; }
  
 for(int i=0; i<lout; i++)
 { d_out[i] = 2.0*pi*(*Vessel)[0].contour(i, 0, "outer")/16.0;
   	if( (*Vessel)[0].contour(i, 0, "outer") > ref_pts_i[2][0] ) d_out[i] = 0.0; }
 	(*Vessel)[0].get_depth(d_in, d_out);
 

// Outboard vacuum vessel
 for(int i=0; i<lin; i++)
 { d_in[i] = 2.0*pi*(*Vessel)[1].contour(i, 0, "inner")/16.0;
   	if( (*Vessel)[1].contour(i, 0, "inner") <= ref_pts_i[2][0] ) d_in[i] = 0.0; }

 for(int i=0; i<lout; i++)
 { d_out[i] = 2.0*pi*(*Vessel)[1].contour(i, 0, "outer")/16.0;
   	if( (*Vessel)[1].contour(i, 0, "outer") <= ref_pts_i[2][0] ) d_out[i] = 0.0; }
 	
 	(*Vessel)[1].get_depth(d_in, d_out);

// Volume and Surface Areas
	for(int i=0; i<2; i++) (*Vessel)[i].get_volume();
	for(int i=0; i<2; i++) (*Vessel)[i].get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof



// ====================================================================================
// Manifolding - SCLL
// Region to the right of OB divertor and above/below OB blankets to route He & PbLi
// into/out of blankets/vessel (LCC copied Repl HT shield to make this section)

void design_manifolding(DesignPoint *Machine, Part (*Manifolding)[2], Part (*RefRHT)[2],
						Part *Ref1, Part *Ref2, Part *Ref3)
						// RefRHT[2] is repl HT shield [0] top
                      	// Ref1 is FW
                      	// Ref2 is blanket I
                      	// Ref3 is blanket II
{ if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour of the Manifolding (not really inner...)
     double ksihi  = (*RefRHT)[0].bounding_box(2, "outer"); // RHS of repl HT shield
     double etahi  = (*RefRHT)[0].bounding_box(3, "outer"); // top of repl HT (y)
     double etahi2 = (*Ref1).bounding_box(3, "outer"); // top of FW (y)
     double ksiFW  = (*Ref1).bounding_box(0, "outer"); // IB, LHS of OB FW (x)
     double ksibl  = (*Ref3).bounding_box(3, "outer"); // Top-right of OB blkt II (y)
//cout << "ksibl = " << ksibl0 << "\n";

// ... ... reference points ... ...
    double ref_pts_i[100][2], center_i[100][2];
    int nptsi=1;
    int resol_i[100];
    string connectors_i[100];

     double dwall = (*Machine).inboard_generic_build(1); // FW thickness

/*     double dsrinb 	   = (*Machine).inboard_generic_build(10); // IB skel ring thickness
     double gapinb 	   = (*Machine).inboard_generic_build(11); // IB skel ring and manifold gap (SR was HT shield)
     double dvinb 	   = (*Machine).inboard_generic_build(12); // IB manifold thickness
     double gapoutb    = (*Machine).outboard_generic_build(16);// OB skel ring and manifold gap (SR was HT shield)
     double dvoutb 	   = (*Machine).outboard_generic_build(17);// OB manifold thickness
     double drephtvert = (*Machine).vertical_generic_build(3); // vert repl HT shield thickness
     double dsrvert    = (*Machine).vertical_generic_build(7); // vert skel ring thickness
     double gapvert    = (*Machine).vertical_generic_build(8); // vert skel ring and manifold gap (SR was HT shield)
     double dvvert 	   = (*Machine).vertical_generic_build(9); // vert manifold thickness
*/

// ... Define the outer contour of entire manifolding... ...
// ... ... reference points ... ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 5;
     double xs[6], ys[6];

	// 1 Start at top of FW
	  xs[0] = ksiFW;
      ys[0] = etahi2;

	// 2 Top-left corner at repl HT shield
      xs[1] = ksiFW;
      ys[1] = etahi;

	// 3 Top-right corner
      xs[2] = ksihi; // 5.7
      ys[2] = etahi;

	// 4 Bottom-right corner
 	  xs[3] = ksihi; // 5.7
      ys[3] = ksibl; // 3.4

     // 5 Top-right of FW
      xs[4] = xs[0]+dwall;
      ys[4] = ys[0];

     // 6 End at start point
      xs[5] = xs[0];
      ys[5] = ys[0];

/*cout << "Manifolding: \n";
cout << "xs[0], ys[0] = " << xs[0] << ", " << ys[0] << "\n";
cout << "xs[1], ys[1] = " << xs[1] << ", " << ys[1] << "\n";
cout << "xs[2], ys[2] = " << xs[2] << ", " << ys[2] << "\n";
cout << "xs[3], ys[3] = " << xs[3] << ", " << ys[3] << "\n";
cout << "xs[4], ys[4] = " << xs[4] << ", " << ys[4] << "\n";
cout << "xs[5], ys[5] = " << xs[5] << ", " << ys[5] << "\n";
*/

	for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    	{ resol_o[i] = 20; } // this is the # pts in between the ref points

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    	connectors_o[i] = "line"; 


// ... Update the Manifolding ...
	(*Manifolding)[0].get_contour_descriptions(nptsi, 	nptso, 	resol_i, 	resol_o,
								ref_pts_i,	ref_pts_o,	connectors_i,	connectors_o);


// Define the Inner Target Curve For Manifolding Shape
	double targcurve[10000][2];  int tgl = 0;

// Connect the Reference Points into Contours
	(*Manifolding)[0].connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
	const double pi = 4.0*atan(1.0);
 	double * d_in; double * d_out;
 	int lin  = (*Manifolding)[0].contour_size("inner"); // 0
 	int lout = (*Manifolding)[0].contour_size("outer"); // 80
//cout << "lout = " << lout << "\n";

 d_in  = new double [lin];
 d_out = new double [lout];

// Loop of the manifolding
 for(int i=0; i<lin;  i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++)  d_out[i] = 2.0*pi*(*Manifolding)[0].contour(i, 0, "outer")/16.0;
 
// Volume and Surface Areas
	(*Manifolding)[0].get_depth(d_in, d_out);
	(*Manifolding)[0].get_volume();
	(*Manifolding)[0].get_surface_areas();
//cout << "Manifolding[0] volume = " << (*Manifolding)[0].volume() << "\n";


//--------------------
// Bottom Manifolding
// ... Define the Contours of Bottom Shield ...
(*Manifolding)[1].get_contour_size(lin,  "inner");
(*Manifolding)[1].get_contour_size(lout, "outer");

 double (*ctr_in)[2], (*ctr_out)[2];
 ctr_in  = new double [lin][2];
 ctr_out = new double [lout][2];

for(int i=0; i<lin; i++)
{ ctr_in[i][0] = 0.0;
  ctr_in[i][1] = 0.0;
  d_in[i] 	   = 0.0; }

for(int i=0; i<lout; i++)
{  ctr_out[i][0] =  (*Manifolding)[0].contour(lout-1-i, 0, "outer");
   ctr_out[i][1] = -(*Manifolding)[0].contour(lout-1-i, 1, "outer"); // mirror across y=0
   d_out[i] 	 = 2.0*pi*ctr_out[i][0]/16.0; 
//cout << ctr_out[i][0] << "	" << ctr_out[i][1] << "\n";
}

// ... Volume and Surface Areas for bottom shield ............
(*Manifolding)[1].get_contours(ctr_in, ctr_out, d_in, d_out);
(*Manifolding)[1].get_volume();
(*Manifolding)[1].get_surface_areas();
//cout << "Manifolding[1] volume = " << (*Manifolding)[1].volume() << "\n";

delete [] d_in;
delete [] d_out;
delete [] ctr_in;
delete [] ctr_out;

} // eof



// ====================================================================================
// 1/4 Vacuum Vessel IB - DCLL

void design_vacuum_vessel_inboard(DesignPoint * Machine, Part (*Vessel)[4], Part * RefBody )
{ 								// RefBody is Manifold with manifold on, else SkeletonRing
if( (*Machine).data_point_filter(0) > 0 ) return; 

// ... Define the Inner Contour ...
     double ksilo  = (*RefBody).bounding_box(0, "outer"); // IB manifold (skel ring)
     double ksihi  = (*RefBody).bounding_box(2, "outer"); // OB manifold (skel ring)
     double etahi  = (*RefBody).bounding_box(3, "outer"); // vert height of manfiold (skel ring)
//cout << "ksilo  = " << ksilo << "\n";
//cout << "ksihi  = " << ksihi << "\n";
//cout << "etahi  = " << etahi << "\n";


// ... ... reference points ... ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];
     
// Define the Outer Contour (DCLL)
	double gapinb, dvinb, dringvert, dsrvert;
if (manifold == "on") {     
       gapinb    = (*Machine).inboard_generic_build(13);  // IB manifold and VV gap
       dsrvert   = (*Machine).vertical_generic_build(7);  // vert skel ring thickness  
       dringvert = (*Machine).vertical_generic_build(9);  // vert manifold thickness  
	 }	
else { gapinb    = (*Machine).inboard_generic_build(13);  // IB skel ring and VV gap
       dringvert = (*Machine).vertical_generic_build(7);  // vert skel ring thickness  
	 }
       dvinb     = (*Machine).inboard_generic_build(14);  // IB VV thickness
	 
	 
/* ALL build data that is used in VV 1-4, 	
	 double dhtinb     = (*Machine).inboard_generic_build(8);   // IB HT shield thickness
     double dringinb   = (*Machine).inboard_generic_build(10);  // IB skel ring thickness
     double gapinb     = (*Machine).inboard_generic_build(13);  // IB skel ring and VV gap
     double dvinb      = (*Machine).inboard_generic_build(14);  // IB VV thickness
     double gapoutb    = (*Machine).outboard_generic_build(18); // OB skel ring and VV gap
     double dvoutb     = (*Machine).outboard_generic_build(19); // OB vacuum vessel thickness
     double dhtrepvert = (*Machine).vertical_generic_build(3);  // vert repl HT shield thickness
     double dringvert  = (*Machine).vertical_generic_build(7);  // vert skel ring thickness  
     double gapvert    = (*Machine).vertical_generic_build(10);	// vert skel ring and VV gap
     double dvvert     = (*Machine).vertical_generic_build(11);	// vert VV thickness
*/
     
// ... ... reference points ... ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 4;
     double xs[10], ys[10];

// Bottom-left on y=0 line
     xs[0] = ksilo - gapinb - dvinb; // ok
     ys[0] = 0.0;

// Top-left (top VV mounts to this)
     xs[1] = xs[0];
     ys[1] = etahi - dringvert - 0.5*dsrvert + gapinb; // ok, 

// Top-right (top VV mounts to this)
     xs[2] = ksilo - gapinb; // ok
     ys[2] = ys[1];

// Bottom-right on y=0 line
     xs[3] = xs[2];
     ys[3] = 0.0;
     
/*cout << "xs[0] = " << xs[0] << "\n";
cout << "ys[0] = " << ys[0] << "\n";
cout << "xs[1] = " << xs[1] << "\n";
cout << "ys[1] = " << ys[1] << "\n";
cout << "xs[2] = " << xs[2] << "\n";
cout << "ys[2] = " << ys[2] << "\n";
cout << "xs[3] = " << xs[3] << "\n";
cout << "ys[3] = " << ys[3] << "\n";
*/

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }


// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-2;

// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line"; }


// Update the Inboard Vacuum Vessel
(*Vessel)[0].get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                      connectors_i, connectors_o);

// Define the Inner Target Curve For Inboard Vacuum Vessel Shape
double targcurve[10000][2]; int tgl=0;

// Connect the Reference Points into Contours
(*Vessel)[0].connect(center_i, center_o, targcurve, tgl);

// Establish Contour Depth for volume calculation
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin  = 1;
 int lout = (*Vessel)[0].contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Vessel)[0].contour(i, 0, "outer")/16.0;

 (*Vessel)[0].get_depth(d_in, d_out);


// Volume and Surface Areas
(*Vessel)[0].get_volume();
(*Vessel)[0].get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// 2/4 Vacuum Vessel OB - DCLL

void design_vacuum_vessel_outboard(DesignPoint * Machine, Part (*Vessel)[4], Part * RefBody )
{ 								// RefBody is Manifold with manifold on, else SkeletonRing
if( (*Machine).data_point_filter(0) > 0 ) return; 

// ... Define the Inner Contour ...
     double ksilo  = (*RefBody).bounding_box(0, "outer"); // IB (X) Skel ring
     double ksihi  = (*RefBody).bounding_box(2, "outer"); // OB (X) Skel ring
     double etahi  = (*RefBody).bounding_box(3, "outer"); // vert height (Y) of Skel ring
     double reflen = (*RefBody).contour_size("outer");    // length of Skel ring contour

// ... ... reference points ... ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];

// Define the Outer Contour (DCLL)
double manoutb, gapoutb, dvoutb, dringvert, gapvert, dvvert;
if (manifold == "on") {
 	   manoutb   = (*Machine).outboard_generic_build(17); // OB manifold thickness
 	   gapoutb   = (*Machine).outboard_generic_build(18); // OB manifold and VV gap
       dringvert = (*Machine).vertical_generic_build(9);  // vert manifold thickness  
       gapvert   = (*Machine).vertical_generic_build(10); // vert manifold and VV gap
	 }
else { gapoutb   = (*Machine).outboard_generic_build(18); // OB skel ring and VV gap
       dringvert = (*Machine).vertical_generic_build(7);  // vert skel ring thickness  
       gapvert   = (*Machine).vertical_generic_build(10); // vert skel ring and VV gap
	 }
       dvoutb    = (*Machine).outboard_generic_build(19); // OB vacuum vessel thickness
	   dvvert    = (*Machine).vertical_generic_build(11); // vert VV thickness
    
     
// ... ... reference points ... ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 4;
     double xs[10], ys[10];

// Midplane of OB VV
     xs[0] = ksihi + gapoutb; // removed + dvoutb;
     ys[0] = 0.0;

// Follow Skeleton ring contour
      for(int i=0; i<reflen; i++)
      { if( (*RefBody).contour(i, 1, "outer") >= (etahi - dvvert) )
       { xs[1] = (*RefBody).contour(i, 0, "outer") - manoutb; } } // ok, removed + dvoutb
         ys[1] = etahi + gapvert; // ok, gap between upper VV and OB VV = 0

// Top-left
     xs[2] = xs[1] + dvoutb; // ok
     ys[2] = ys[1];

// Top-right
     xs[3] = xs[0] + dvoutb; // ok
     ys[3] = 0.0;

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }


// ... Close the loop by adding a symmetric profile in the clockwise direction:
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-2;

// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "target";
      if(i==1 || i==4) connectors_o[i] = "line"; }

// Update the Outboard Vacuum Vessel
(*Vessel)[1].get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                      connectors_i, connectors_o);


// Define the Inner Target Curve For Outboard Vacuum Vessel Shape
double targcurve[10000][2]; 
int tgl = (*RefBody).contour_size("outer");
  for(int i=0; i<tgl; i++)
  {   targcurve[i][0] = (*RefBody).contour(i, 0, "outer");
      targcurve[i][1] = (*RefBody).contour(i, 1, "outer"); }


// Connect the Reference Points into Contours
(*Vessel)[1].connect(center_i, center_o, targcurve, tgl);

// Establish Contour Depth for volume calculation
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin  = 1;
 int lout = (*Vessel)[1].contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Vessel)[1].contour(i, 0, "outer")/16.0;

 (*Vessel)[1].get_depth(d_in, d_out);


// Volume and Surface Areas

(*Vessel)[1].get_volume();
(*Vessel)[1].get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// 3/4 Vacuum Vessel Bottom DCLL - "By Symmetry"

void copy_vacuum_vessel_top_into_bottom_by_symmetry(Part (*Vessel)[4])
{
  int lin  = (*Vessel)[2].contour_size("inner");
  int lout = (*Vessel)[2].contour_size("outer");

 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;

  d_in  = new double [lin];
  d_out = new double [lout];

(*Vessel)[3].get_contour_size(lin, "inner");
(*Vessel)[3].get_contour_size(lout, "outer");

double (*ctr_in)[2], (*ctr_out)[2];

ctr_in  = new double [lin][2];
ctr_out = new double [lout][2];

for(int i=0; i<lin; i++)
{ ctr_in[i][0] = 0.0;
  ctr_in[i][1] = 0.0;
  d_in[i]      = 0.0; }

for(int i=0; i<lout; i++)
{  ctr_out[i][0] =  (*Vessel)[2].contour(lout-1-i, 0, "outer");
   ctr_out[i][1] = -(*Vessel)[2].contour(lout-1-i, 1, "outer");
   d_out[i]      = 2.0*pi*ctr_out[i][0]/16.0; }

(*Vessel)[3].get_contours(ctr_in, ctr_out, d_in, d_out);


// ... Volume and Surface Areas ...
(*Vessel)[3].get_volume();
(*Vessel)[3].get_surface_areas();

delete [] d_in;
delete [] d_out;
delete [] ctr_in;
delete [] ctr_out;

} // eof


// ====================================================================================
// 4/4 Vacuum Vessel Top DCLL - "By Design"

void design_vacuum_vessel_top(DesignPoint * Machine, Part (*Vessel)[4], Part * RefBody )
{ 								// RefBody is Manifold with manifold on, else SkeletonRing
if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour ...
     double ksilo  = (*RefBody).bounding_box(0, "outer"); // IB of manifold (SR)
     double ksihi  = (*RefBody).bounding_box(2, "outer"); // OB of manifold (SR)
     double etahi  = (*RefBody).bounding_box(3, "outer"); // top of manifold (SR)
     double reflen = (*RefBody).contour_size("outer");
//cout << "ksilo = " << ksilo << "\n";
//cout << "ksihi = " << ksihi << "\n";
//cout << "etahi = " << etahi << "\n";

// ... ... reference points ... ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];

// Define the Outer Contour (DCLL)
	double dhtinb, dringinb, gapinb, dvinb, dvoutb, dhtrepvert, dringvert, 
		gapvert, dvvert, dsrvert;
if (manifold == "on") {
	   dhtinb     = (*Machine).inboard_generic_build(8);   // IB HT shield thickness = 0
       dringinb   = (*Machine).inboard_generic_build(12);  // IB manifold thickness
       gapinb     = (*Machine).inboard_generic_build(13);  // IB manifold and VV gap
       dringvert  = (*Machine).vertical_generic_build(9);  // vert manifold thickness  
       gapvert    = (*Machine).vertical_generic_build(10); // vert manifold and VV gap
	 }
else { dhtinb     = (*Machine).inboard_generic_build(8);   // IB HT shield thickness
       dringinb   = (*Machine).inboard_generic_build(10);  // IB skel ring thickness
       gapinb     = (*Machine).inboard_generic_build(13);  // IB skel ring and VV gap
       dringvert  = (*Machine).vertical_generic_build(7);  // vert skel ring thickness  
       gapvert    = (*Machine).vertical_generic_build(10); // vert skel ring and VV gap
	 }
       dvinb      = (*Machine).inboard_generic_build(14);  // IB VV thickness
       dvoutb     = (*Machine).outboard_generic_build(19); // OB vacuum vessel thickness
       dhtrepvert = (*Machine).vertical_generic_build(3);  // vert repl HT shield thickness
       dsrvert    = (*Machine).vertical_generic_build(7);  // vert skel ring thickness  
       dvvert     = (*Machine).vertical_generic_build(11); // vert VV thickness
	
	
// ... ... reference points ... ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 8;
     double xs[10], ys[10];

     xs[0] = ksilo - gapinb - dvinb; // top-left on top of IB VV, then go clock-wise
     ys[0] = etahi - 0.5*dsrvert - dringvert + gapvert; // removed - 0.5*dringvert
// Should really be referencing the top of IB VV rather than this +/- game...

     xs[1] = ksilo + dhtrepvert + dhtinb + dringinb + 0.4; // top-left of inner contour, temp
     ys[1] = etahi + dvvert; // removed + gapvert

// Connects top VV to OB top-right
     for(int i=0; i<reflen; i++) 
      { if( (*RefBody).contour(i, 1, "outer") >= (etahi - dvvert) ) 
       { xs[2] = (*RefBody).contour(i, 0, "outer") - gapvert; } } // top-right, removed + dvoutb, ok
         ys[2] = etahi + gapvert + dvvert; 

     xs[3] = xs[2] + dvoutb; // bottom-right
     ys[3] = ys[2] - dvvert;

     xs[4] = xs[2]; // intersection w/ top-left of OB VV
     ys[4] = ys[3];

     xs[5] = xs[1] + gapvert; // bottom-left
     ys[5] = ys[3];

     xs[6] = xs[0] + dvinb; // bottom-left of slope
     ys[6] = ys[0];

     xs[7] = xs[0] + 1.e-9; // back to beginning
     ys[7] = ys[0] + 1.e-9;

/*cout << "xs[0], ys[0] = " << xs[0] << ", " << ys[0] << "\n";
cout << "xs[1], ys[1] = " << xs[1] << ", " << ys[1] << "\n";
cout << "xs[2], ys[2] = " << xs[2] << ", " << ys[2] << "\n";
cout << "xs[3], ys[3] = " << xs[3] << ", " << ys[3] << "\n";
cout << "xs[4], ys[4] = " << xs[4] << ", " << ys[4] << "\n";
cout << "xs[5], ys[5] = " << xs[5] << ", " << ys[5] << "\n";
cout << "xs[6], ys[6] = " << xs[6] << ", " << ys[6] << "\n";
cout << "xs[7], ys[7] = " << xs[7] << ", " << ys[7] << "\n";
*/

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 100; }

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
      if(i==0 || i==5) connectors_o[i] = "target"; }

// Update the Inboard Vacuum Vessel
(*Vessel)[2].get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, ref_pts_o,
                                 	  connectors_i, connectors_o);

// Define the Inner Target Curve For Outboard Vacuum Vessel Shape
double targcurve[10000][2]; 
int tgl = (*RefBody).contour_size("outer");
  for(int i=0; i<tgl; i++)
  {   targcurve[i][0] = (*RefBody).contour(i, 0, "outer");
      targcurve[i][1] = (*RefBody).contour(i, 1, "outer"); }

// Connect the Reference Points into Contours
(*Vessel)[2].connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation
 const double pi = 4.0*atan(1.0); // pi = 3.14159....
 double * d_in; double * d_out;
 int lin  = 1;
 int lout = (*Vessel)[2].contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Vessel)[2].contour(i, 0, "outer")/16.0;

 (*Vessel)[2].get_depth(d_in, d_out);


// Volume and Surface Areas
(*Vessel)[2].get_volume();
(*Vessel)[2].get_surface_areas();

delete [] d_in;
delete [] d_out;

} // eof



// START OF COIL DESIGNS
// ====================================================================================
// ====================================================================================
// TF Coil

// ... r-theta plane, "PPPL Version, Conservative Costing"

void design_tf_coil_rtheta_PPPLConservative(DesignPoint * Machine, Part * Magnet, Part * Plasma, double Bt)
{
	const double pi  = 4.0*atan(1.0);
	const double muo = 4.0*pi*1.e-7;

/* c
c...toroidal field coil
c
c...ntfcoil    =  number of TF coils
c...btmaxlim   =  maximum allowed toroidal field at the TF magnet
c...xjsclim    =  maximum allowed superconducting current density in TF
c...fstarea    =  TF structural material area fraction
c...fcoolarea  =  TF coolant area fraction
c...finsarea   =  TF insulation area fraction
c...fcuarea    =  TF copper stabilizer area fraction
c...sigall     =  stress allowable for TF coil
c...sigcu      =  stress allowable in Cu in TF coil
c...xjcu       =  maximum current density in Cu stabilizer
c...rratio     =  major radius ratio (midplane outboard over midplane inboard)
c...sigtf      =  total TF coil stress allowable
c...tfcur      =  current in each TF coil
c...rtfout     =  outer major radius of inboard TF coil
c...btmax      =  maximum toroidal field at TF coil
c...fscarea    =  TF superconducting area fraction
c...tfarea     =  area of annulus of inboard TF coil legs
c...rtfin      =  inner major radius of inboard TF coil  */

// References for below???
      const double ntfcoil = 16.0; // number of TF coils
//      const double btmaxlim = 21.0; // Bt is allowed to go big - it is post-filtered
      const double xjsclim = 2.5e8;
      const double fstarea = 0.65;
      const double fcoolarea = 0.03;
      const double finsarea = 0.07;
      const double fcuarea = 0.15;
//      const double sigall = 8.e8;
//      const double sigcu = 8.e8;
//      const double xjcu = 2.e8;
//      const double rratio = 3.5;
      const double sigtf = 5.e8;
      double fscarea  = 1. - fstarea - fcuarea - fcoolarea - finsarea;
      double tfcur    = Bt * 2.0 * pi * (*Machine).plasma_major_radius() /
      					(1.2566e-06 * ntfcoil);   // current in each TF coil
      double deltainb = 0.0;

// Sum of IB radial builds from plasma surface until TF coil
// LCC WATCH this scan in case radial build grows!
      for( int i=0; i<9; i++)
      { deltainb = deltainb + (*Machine).inboard_generic_build(i); }

      double rtfout = (*Machine).plasma_major_radius() -
                      (*Machine).plasma_minor_radius() -
                      deltainb; // outer major radius of inboard tf coil

      double btmax = 1.2566e-06 * ntfcoil * tfcur / (2.0 * pi * rtfout);

      double rtfin = pow( pow(rtfout, 2.0) -
                          (  ( tfcur/(pi*xjsclim) )*( (1-fstarea)/(fscarea) )  ) -
                          ( pow(rtfout, 2.0)*pow(btmax, 2.0) )/(sigtf*1.2566e-6),
                          0.5);

      double dtfi = rtfout - rtfin;



// Get TF Cross-Section Material Components
     double delta = dtfi;
     double width = 2.0*pi*rtfout/ntfcoil;
     double Itf = 40000.0;  // TF coil current

     (*Magnet).get_cross_section(0, delta, width);
     (*Magnet).get_material_component("SS316", fstarea, 1.0, 0); // was Structconserv
     
     // Number of turns per winding pack
     int Nwp = int( Bt*2.0*pi*(*Machine).plasma_major_radius()/(ntfcoil*Itf*muo) ); 
	 if(Nwp ==0) Nwp = 1;
(*Machine).get_number_of_tf_cable_turns(Nwp);
double Rinnermost = rtfin;
double Routermost = 1.03*3.87*rtfout;
(*Machine).get_tf_coil_radii(Rinnermost, Routermost); // gets innermost and outermost radius of the TF coil

} // eof


// ====================================================================================
// ... r-theta plane, "PPPL Version"

void design_tf_coil_rtheta_PPPL(DesignPoint * Machine, Part * Magnet, Part * Plasma, double Bt)
{
	const double pi = 4.0*atan(1.0);
	const double muo = 4.0*pi*1.e-7;

/* The following are from C. Kessel's code "syseng.f"
c...toroidal field coil
c
c...ntfcoil  	=  number of TF coils
c...btmaxlim  	=  maximum allowed toroidal field at the TF magnet
c...xjsclim  	=  maximum allowed superconducting current density in TF
c...fstarea  	=  TF structural material area fraction
c...fcoolarea  	=  TF coolant area fraction
c...finsarea  	=  TF insulation area fraction
c...fcuarea  	=  TF copper stabilizer area fraction
c...sigall  	=  stress allowable for TF coil
c...sigcu  		=  stress allowable in Cu in TF coil
c...xjcu  		=  maximum current density in Cu stabilizer
c...rratio  	=  major radius ratio (midplane outboard over midplane inboard)
c...sigtf  		=  total TF coil stress allowable
c...tfcur  		=  current in each TF coil
c...rtfout  	=  outer major radius of inboard TF coil
c...btmax  		=  maximum toroidal field at TF coil
c...fscarea  	=  TF superconducting area fraction
c...tfarea  	=  area of annulus of inboard TF coil legs
c...rtfin  		=  inner major radius of inboard TF coil  */

      const double ntfcoil = 16.0; // number of TF coils
//      const double btmaxlim = 21.0;
      const double xjsclim = 2.5e8;
      const double fstarea = 0.65;
      const double fcoolarea = 0.03;
      const double finsarea = 0.07;
      const double fcuarea = 0.15;
//      const double sigall = 8.e8;
//      const double sigcu = 8.e8;
//      const double xjcu = 2.e8;
//      const double rratio = 3.5;
      const double sigtf = 5.e8;
      
      double fscarea = 1. - fstarea - fcuarea - fcoolarea - finsarea;
      double tfcur = Bt*2.0*pi*(*Machine).plasma_major_radius()/(1.2566e-06*ntfcoil);   // current in each tf coil      
      double deltainb = 0.0;

      for( int i=0; i<9; i++)
      { deltainb = deltainb + (*Machine).inboard_generic_build(i); }
      // sum of inboard radial builds from plasma surface until tf coil

      double rtfout = (*Machine).plasma_major_radius() -
                      (*Machine).plasma_minor_radius() -
                      deltainb; // outer major radius of inboard tf coil
      
      double btmax = 1.2566e-06*ntfcoil*tfcur/(2.*pi*rtfout);
      
      double rtfin = pow(pow(rtfout, 2.0) -
                     ((tfcur/(pi*xjsclim))*((1-fstarea)/(fscarea))) -
                     (pow(rtfout, 2.0)*pow(btmax, 2.0) )/(sigtf*1.2566e-6), 0.5);
      double dtfi = rtfout - rtfin;
      


// Get TF Cross-Section Material Components
     double delta = dtfi;
     double width = 2.0*pi*rtfout/ntfcoil;
     double Itf = 40000.0;

   (*Magnet).get_cross_section(0, delta, width);
   (*Magnet).get_material_component("JK2LB", fstarea, 1.0, 0);
   (*Magnet).get_material_component("Nb3Sn", fscarea, 1.0, 1);

   int Nwp = int( Bt*2.0*pi*(*Machine).plasma_major_radius()/(ntfcoil*Itf*muo) ); // number of turns per winding pack
      
   if(Nwp ==0) Nwp = 1;
(*Machine).get_number_of_tf_cable_turns(Nwp);
double Rinnermost = rtfin;
double Routermost = 1.03*3.87*rtfout;
(*Machine).get_tf_coil_radii(Rinnermost, Routermost); // gets innermost and outermost radius of the TF coil

} // eof




// ====================================================================================
// ********* CURRENT TF r-theta algorithm IN USE **************************************
// ====================================================================================
// ... r-theta plane, "Next Step", conservative costing (see "TF Magnet.doc" in ASC Docs)
// This algorithm is selected in Aries.cpp ~423 "TF Coil"

// TF Coil - both, (DCLL is same except for (*VVRefBody)[4] vs [2]
void design_tf_coil_rtheta_NextStepConservative(DesignPoint *Machine,
		Part *Magnet, Part *Plasma, Part (*VVRefBody)[4], double Bt)
{ 
if((*Machine).data_point_filter(0) > 0) return;
   
   const double pi  = 4.0*atan(1.0);
   const double muo = 4.0*pi*1.0e-7;
   
// Open the input file to read the magnet parameters:
  ifstream fin_mag;  fin_mag.open("InputData/magnets.data");
  fin_mag.width(10); fin_mag.precision(2);
  fin_mag.setf(ios_base::uppercase, ios_base::scientific);
  int itxt_mag = 57; // length of variable description in the file
  int inum_mag = 13; // length of the formated number in the file
  int iskip_mag = itxt_mag + inum_mag; // number of places occupied per variable


// ... maximum stress in protective sheath around strands
   double sigma_sheath; fin_mag.seekg(0*iskip_mag + itxt_mag + 1); fin_mag >> sigma_sheath;

// ... Maximum allowed stress in the coil casing
   double sigmab; fin_mag.seekg(1*iskip_mag + itxt_mag + 1); fin_mag >> sigmab;

// Define the bounding box for the plasma
   double ksiloP = (*Plasma).bounding_box(0, "outer");
   double ksihiP = (*Plasma).bounding_box(2, "outer");
//cout << "ksiloP   = " << ksiloP << "\n";
//cout << "ksihiP   = " << ksihiP << "\n";

// Define the bounding box for the Vacuum Vessel
   double ksihiVV = (*VVRefBody)[1].bounding_box(2, "outer");
   double ksiloVV = (*VVRefBody)[0].bounding_box(0, "outer");
//cout << "ksihiVV  = " << ksihiVV << "\n";
//cout << "ksiloVV  = " << ksiloVV << "\n";

// Define the IB radius of the TF coil that corresponds to its inner surface. 
// The IB thicknesses of all the power core elements enclosed within the TF coil 
// including the gaps in between are subtracted from the ksiloP of the plasma.

// Gap between the VV and IB TF coil
	double gapvvtf  = (*Machine).inboard_generic_build(17); // VV and IB TF coil gap
	double rTFinner = ksiloVV - gapvvtf; // radius of inner TF coil
//cout << "gapvvtf  = " << gapvvtf << "\n";
//cout << "rTFinner = " << rTFinner << "\n";

// Maximum Magnetic Field
   double Ro 	= (*Machine).plasma_major_radius(); // plasma_major_radius
   double Btmax = Bt*Ro/rTFinner;
   (*Machine).get_maximum_magnetic_field_at_tf_coil(Btmax); // puts BTmax into this fcn
//cout << "Bt = " << Bt << "\n";

// Super-conducting current table from ZD TF_Magnet.doc
// Current Density in Superconductor versus Btmax function for Nb3Sn @ 4.2 K
   double btmaxVec[5] = {6.0,	 10.0, 	 14.0, 	 16.0,   18.0  };
   double jscVec[5] =   {1.52e9, 1.35e9, 8.84e8, 5.43e8, 1.30e8};

   double Aturn, Astrands, Acu, Asc, Asheath, Ains;
   double delta_inner, delta_outer, delta_top, wTF, dTF, b;
   int Nwp;

   if(Btmax < btmaxVec[0] || Btmax > btmaxVec[4])
   { cout << "Btmax (" << Btmax << ") out of range\n";  // Out of range of table data
     (*Machine).get_data_point_filter(1, 1); return; }

   double dist = fabs(Btmax - btmaxVec[0]);
   int imark = 0;
   for(int i=0; i<5; i++)
   { if( fabs(Btmax - btmaxVec[i]) < dist)
   	{ dist = fabs(Btmax - btmaxVec[i]);  // This just finds the nearest entry from the table
     imark = i; } }
//   double Jsc = jscVec[imark]; // nearest current density value of the superconductor


// Super-conducting current polynomial
	// LCC- Use polynomial equation INSTEAD of above table, based on poly fit to table values
	// C. Kessel also uses this in his code "sysengr.f" ~364
	double Jsc = -9.1226*pow(Btmax,2)*1.0e6 + 1.0301*Btmax*1.0e8 + 1.2308*1.0e9;
	

// Cross section of the winding pack:
// ... Current density in copper
	double Jcu; fin_mag.seekg(2*iskip_mag + itxt_mag + 1); fin_mag >> Jcu; // 2.4e8
   
   
// Reduce Magnet Current by 1.5, 3
	// Current estimates may be too optimistic 
//	Jsc = Jsc/3.0;
//	Jcu = Jcu/3.0;
//cout << "Btmax = " << Btmax << "\n";
//cout << "Jsc = " << Jsc << "\n";
//cout << "Jcu = " << Jcu << "\n";


// ... Volume fraction of He coolant
   double f_he; fin_mag.seekg(3*iskip_mag + itxt_mag + 1); fin_mag >> f_he; // 25%

// ... Sheath thickness
   double dsheath = pow(Btmax,2)/(2.0*muo*sigma_sheath); // OK

// ... Sheath fraction
   double f_sheath = 1.0 - pow(1.0 - dsheath, 2); // OK

// ... Volume fraction of the insulator
   double f_ins; fin_mag.seekg(4*iskip_mag + itxt_mag + 1); fin_mag >> f_ins; // 10%

// ... Nominal current in TF coil
   double Itf; fin_mag.seekg(5*iskip_mag + itxt_mag + 1); fin_mag >> Itf; // 40,000 A
   
   
// Cross-sectional areas:
	Acu      = Itf/Jcu;
	Asc      = Itf/Jsc;
	Astrands = Acu + Asc;
	Aturn    = Astrands/(1.0 - f_he - f_sheath - f_ins); // C.K. sysengr.f line ~376
	// Cross-sectional area of a single turn in a winding pack
	Ains     = Aturn*f_ins;
	Asheath  = Aturn*f_sheath;


// Coil Casing Cross Section

// Number of TF coils, now a variable in magnets.data input file, NEEDS more widespread use
   int Ntf; fin_mag.seekg(14*iskip_mag + itxt_mag + 1); fin_mag >> Ntf; // 16 TF coils
//cout << Ntf << "\n";

// ... Side thickness of the coil casing is b/2
   fin_mag.seekg(6*iskip_mag + itxt_mag + 1); fin_mag >> b; // b = 0.02

// Area of a TF coil single winding pack
   //double ATF = Aturn*pow(Bt,2.0)*pi*Ro/(Ntf*Itf*muo);
   double ATF = Aturn*Bt*2.0*pi*Ro/(Ntf*Itf*muo); // OK

// Number of turns per winding pack (WP), ~236   
   Nwp = int(ATF/Aturn); 
   (*Machine).get_number_of_tf_cable_turns(Nwp); // puts Nwp into this fcn

//cout << "ATF = " << ATF << "\n";
//cout << "Aturn = " << Aturn << "\n";


// Initial guess for casing thicknesses
   	double a, R1, R2, coeff;
	a=b; // b = side thickness of coil casing
for(int j=0; j<5; j++)
{ // start LOOP

// Calculate R1 = inner radius of inner TF coil leg
	R1 = rTFinner - 2.0*a;
	for(int i=0; i<5; i++)
		{ R1 = rTFinner - a - ATF/( 2.0*pi*R1/Ntf - b); } // similar to CK sysengr.f ~382
	wTF = 2.0*pi*R1/Ntf - b;  // width of the TF coil
	dTF = ATF/wTF;  		  // winding pack (conductors) thickness if you squished them all together
	// This is smaller than ITER
	// Alternatives: dTF = rTFinner - a - R1; or = ATF/( 2.0*pi*R1/Ntf - b);

// Calculate R2 based on the clearance on the outboard side of the TF coils
// needed for removal of the power core module for maintenance purpose
// R2 = inner radius of outer TF coil leg

// ... Clearance between the Vacuum Vessel and the TF coil outmost edge
   double cTFVV; fin_mag.seekg(7*iskip_mag + itxt_mag + 1); fin_mag >> cTFVV; // = 0.1 m
   R2 = ksihiVV + R1 + cTFVV/sin(pi/Ntf); 


// Structure calculation based on bending:
   double G1 = pow(3.0/(16.0*sigmab*muo), 0.5);
   double G2 = Bt*Ro/(R2-R1);
   double G3 = pow((R2*R2-R1*R1)*(1+log(R2/R1)) + 
   				    R1*R2*(2-(log(R2/R1))*(log(R2/R1))),0.5 ); // extra log is for log^2
//cout << "G1 = " << G1 << "\n";
//cout << "G2 = " << G2 << "\n";
//cout << "G3 = " << G3 << "\n";

// Coil casing thickness, including structure (casing and internal plates) from "TF Magnet.doc"
	coeff = 0.208; // OK for ACT-I, II
	a = coeff*G1*G2*G3;  
	// 0.208 coeff is a casing thickness multiplier based on scaling FEA results to match AT casing.
	// This is also in CK's syseng.f but ASC Docs 'TF magnet.doc' has .3351...
	// Is a the total coil casing, so a/2 + winding pack + a/2 = total TF thickness? Yes.
	// This includes coil case and/or all internal structure.
    // There is a tradeoff between internal and coil casing structure Ð if one goes down, the other must go up.

} // END LOOP

//cout << "coefficient = " << coeff << "\n";
//cout << "structure (internal and coil case) thickness, a = " << a << "\n";
//cout << "winding pack only thickness, dTF = " << dTF << "\n";
//cout << "TF coil width, wTF = " << wTF << "\n";
//cout << "total magnet = dTF + a = " << dTF+a << "\n";
//cout << "R1 = " << R1 << "\n";
//cout << "R2 = " << R2 << "\n";

// Defining parameters of the cross section
// (the thickness is uniform along the TF coil perimeter for now)
delta_inner = a + dTF; // total thickness of inner leg equal to casing thickness + magnet thickness
delta_outer = a + dTF; // total thickness of outer leg equal to casing thickness + magnet thickness
delta_top   = a + dTF; // total thickness of top/bottom part equal to casing thickness + magnet thickness
//cout << "delta_inner = " << delta_inner << "\n"; // all 0.3265 m since thickness is uniform


// Update the outer radius of the outer surface of the TF coil:
	R2 = R2 + delta_outer;

// Get the inner-most and outer-most radius of the TF Coil:
   (*Machine).get_tf_coil_radii(R1, R2);
   double delta = max(delta_inner, max(delta_outer, delta_top));
   double width = wTF + b; // width of TF coil

// Cross Sectional Area of TF Coil

// Cross sectional area of the composition
   double Acomp 			= delta*width; 			 // cross sectional area of the composition
   double Astructure 		= delta*width - dTF*wTF; // steel (should be ~75% of WP+CC)
   double vfrac_structure 	= Astructure/Acomp; 	 // volume fraction of SS316 structure
   double vfrac_sc 			= Asc*Nwp/Acomp; 		 // superconductor (Nb3Sn)
   double vfrac_cu 			= Acu*Nwp/Acomp; 		 // copper (Cu)
   double vfrac_sheath 		= Asheath*Nwp/Acomp; 	 // sheath (Inconel)
   double vfrac_ins 		= Ains*Nwp/Acomp; 		 // insulator (Polyamide)
   double vfrac_he = 1.0-vfrac_structure-vfrac_sc-vfrac_cu-vfrac_sheath-vfrac_ins; // He coolant
//cout << "vfrac_structure = " << vfrac_structure << "\n";
//cout << "vfrac_sc = " << vfrac_sc << "\n";
//cout << "vfrac_cu = " << vfrac_cu << "\n";
//cout << "vfrac_sheath = " << vfrac_sheath << "\n";
//cout << "vfrac_ins = " << vfrac_ins << "\n";
//cout << "He fraction = " << 1-vfrac_structure-vfrac_sc-vfrac_cu-vfrac_sheath-vfrac_ins << "\n";


// Get TF Cross-Section Material Components
// Gets the material costs, mass, etc from Part.cpp based on:
// (1) material name, (2) volume fraction, (3) density fraction, and (4) the cost base.
   (*Magnet).get_cross_section(0, delta, width);
   (*Magnet).get_material_component("SS316", 		vfrac_structure, 1.0, 0); // was "Structconserv"
   (*Magnet).get_material_component("Nb3Sn", 		vfrac_sc, 		 1.0, 1);
   (*Magnet).get_material_component("Cu", 			vfrac_cu, 		 1.0, 2);
   (*Magnet).get_material_component("Inconel", 		vfrac_sheath, 	 1.0, 3);
   (*Magnet).get_material_component("Polyamide", 	vfrac_ins, 		 1.0, 4);
   (*Magnet).get_material_component("He",		 	vfrac_he, 		 1.0, 5); // LCC added
   
fin_mag.close();
} // eof

// ====================================================================================
// ======= END of CURRENT TF algorithm ========+=======================================
// ====================================================================================


// ... r-theta plane, "Next Step"

void design_tf_coil_rtheta_NextStep(DesignPoint * Machine, Part * Magnet, Part * Plasma, double Bt)
{
if( (*Machine).data_point_filter(0) > 0 ) return;
   const double pi = 4.0*atan(1.0);
   const double muo = 4.0*pi*1.e-7;
   const double sigma_sheath = 8.0e8;
   const double sigmab = 400.0e6;


// Define the bounding box for the plasma
     double ksiloP = (*Plasma).bounding_box(0, "outer");

// Define the IB radius of the TF coil that corresponds to its inner surface. 
// The IB thicknesses of all the power core elements enclosed within the TF coil
// including the gaps in between are subtracted from the ksiloP of the plasma.

double deltainb = 0.0;
for( int i=0; i<9; i++)
{ deltainb = deltainb + (*Machine).inboard_generic_build(i); }
double rTFinner = ksiloP -  deltainb;


// Maximum Magnetic Field
double Btmax = Bt*(*Machine).plasma_major_radius()/rTFinner;
// cout << "peak field at the coil = " << Btmax << "\n";

// Current Density in Superconductor versus Btmax function for Nb3Sn @ 4.2 K
double btmaxVec[5] = {6.0, 10.0, 14.0, 16.0, 18.0};
double jscVec[5] = {1.52e9, 1.35e9, 8.84e8, 5.43e8, 1.30e8};

// LCC modified jscVec 6/29/10 to determine effect of reducing magnetic current by 1.5, 3
// current estimates may be too optimistic 
//double jscVec[5] = {(1.52e9)/1.5, (1.35e9)/1.5, (8.84e8)/1.5, (5.43e8)/1.5, (1.30e8)/1.5};

double Itf, Aturn, Astrands, Acu, Asc, Asheath, Ains;
double delta_inner, delta_outer, delta_top;
double wTF, dTF, b;
int Nwp;

for(int kbtmax=0; kbtmax<3; kbtmax++)  // Iterative loop over Btmax
{
if(Btmax < btmaxVec[0] || Btmax > btmaxVec[4])
{
  cout << "Btmax out of range\n";
  (*Machine).get_data_point_filter(1, 1);
  return;
}

double dist = fabs(Btmax - btmaxVec[0]);
int imark = 0;

for(int i=0; i<5; i++)
{
if( fabs(Btmax - btmaxVec[i]) < dist)
{
  dist = fabs(Btmax - btmaxVec[i]);
  imark = i;
}
}


// Cross section of the winding pack:
double Jsc = jscVec[imark]; // current density of the superconductor
double Jcu = 2.24e8; 		// current density in copper

// LCC modified Jcu 6/29/10 to determine effect of reducing magnetic current by 1.5, 3
//double Jcu = (2.24e8)/1.5; // current density in copper

double f_he = 0.25; // helium fraction
double dsheath = pow(Btmax,2)/(2.0*muo*sigma_sheath);   // sheath thickness
double f_sheath = 1.0 - pow(1.0-dsheath, 2); // sheath fraction
double f_ins = 0.1; // insulator fraction


// Cross-sectional areas:
Itf = 40000.0;  //current [A]
Acu = Itf/Jcu;
Asc = Itf/Jsc;
Astrands = Acu + Asc;
Aturn = Astrands/(1.0 - f_he - f_sheath - f_ins);
Ains = Aturn*f_ins;
Asheath = Aturn*f_sheath;


// Coil Casing Cross Section
double R2 = 1.03*3.87*rTFinner; // outer external structure radius
int Ntf = 16;  // number of TF coils
b = 0.02; // side thickness of the casing is b/2
double ATF = Aturn*Bt*2.0*pi*(*Machine).plasma_major_radius()/(Ntf*Itf*muo);
Nwp = int(ATF/Aturn); // number of turns per winding pack

(*Machine).get_number_of_tf_cable_turns(Nwp);


// Initial guess for casing thicknesses
double a, R1;
  a=b;
for(int j=0; j<5; j++)
{

// Calculate R1
R1= rTFinner - 2*a;
for(int i=0; i<5; i++)
{
    R1 = rTFinner - a - ATF/( 2*pi*R1/Ntf - b);
}
wTF = 2*pi*R1/Ntf - b;   // width of the TF coil
dTF = rTFinner - a - R1;  // thickness of the TF coil


// Structure Calculation based on bending:
double G1 = pow(3.0/(16.0*sigmab*muo), 0.5);
double G2 = Bt*(*Machine).plasma_major_radius()/(R2-R1);
double G3 = pow( (R2*R2-R1*R1)*(1.0+log(R2/R1)) + R1*R2*(2.0 - (log(R2/R1))*(log(R2/R1))  ),0.5 );
a = 0.208*G1*G2*G3;
}


// Defining parameters of the cross section
delta_inner = a + dTF;
delta_outer = a + dTF;
delta_top   = a + dTF;

Btmax = Bt*(*Machine).plasma_major_radius()/(rTFinner - delta_inner/2.0);


// Get the inner-most and outer-most radius of the TF Coil:
  if(kbtmax == 2) (*Machine).get_tf_coil_radii(R1, R2);
  
} // End of iterative loop over Btmax


double delta = max(delta_inner, max(delta_outer, delta_top));
double width = wTF + b;


// Cross Sectional Area of TF Coil and Volume Fractions

double Acomp = delta*width;    // cross sectional area of the composition
// ... steel
double Astructure = delta*width - dTF*wTF; // cross sectional area of the structure
double vfrac_structure =  Astructure/Acomp;
// ... superconductor
double vfrac_sc = Asc*Nwp/Acomp;
// ... copper
// double vfrac_cu = Acu*Nwp/Acomp;
// ... sheath
// double vfrac_sheath = Asheath*Nwp/Acomp;
// ... insulator
// double vfrac_ins = Ains*Nwp/Acomp;


// Get TF Cross-Section Material Components
   (*Magnet).get_cross_section(0, delta, width);
   (*Magnet).get_material_component("JK2LB", vfrac_structure, 1.0, 0);
   (*Magnet).get_material_component("Nb3Sn", vfrac_sc, 1.0, 1);

} // eof


// ====================================================================================
// ... r-theta plane, HTS

void design_tf_coil_rtheta_HTS(DesignPoint * Machine, Part * Magnet, Part * Plasma, double Bt)
{
if( (*Machine).data_point_filter(0) > 0 ) return;
   const double pi = 4.0*atan(1.0);
   const double muo = 4.0*pi*1.e-7;
   const double sigma_sheath = 1.6e9;
   const double sigmab = 400.0e6;
   

// Define the bounding box for the plasma
     double ksiloP = (*Plasma).bounding_box(0, "outer");
     
// Define the Inboard radius of the TF coil that corresponds to its inner
// surface. The inboard thicknesses of all the power core elements enclosed
// within the TF coil including the gaps in between are subtracted from the
// ksiloP of the plasma.

double deltainb = 0.0;
for( int i=0; i<9; i++)
{
     deltainb = deltainb + (*Machine).inboard_generic_build(i);
}
double rTFinner = ksiloP -  deltainb;


// Maximum Magnetic Field

double Btmax = Bt*(*Machine).plasma_major_radius()/rTFinner;

// cout << "peak field at the coil = " << Btmax << "\n";


// Current Density in Superconductor versus Btmax function for YBCO @ 75K

double btmaxVec[5] = {6.0, 10.0, 14.0, 16.0, 18.0};
double jscVec[5] = {9.00E+08, 6.75E+08, 5.40E+08, 4.86E+08, 4.62E+08};
double Itf, Aturn, Astrands, Acu, Asc, Asheath, Ains;
double delta_inner, delta_outer, delta_top;
double wTF, dTF, b;
int Nwp;

for(int kbtmax=0; kbtmax<3; kbtmax++)  // Iterative loop over Btmax
{
if(Btmax < btmaxVec[0] || Btmax > btmaxVec[4])
{
  cout << "Btmax out of range\n";
  (*Machine).get_data_point_filter(1, 1);
  return;
}
double dist = fabs(Btmax - btmaxVec[0]);
int imark = 0;

for(int i=0; i<5; i++)
{
if( fabs(Btmax - btmaxVec[i]) < dist)
{
  dist = fabs(Btmax - btmaxVec[i]);
  imark = i;
}
}


// Cross section of the winding pack:

double Jsc = jscVec[imark]; // current density of the superconductor
double Jcu = Jsc/2.0; // current density in copper

// LCC modified Jcu 6/29/10 to determine effect of reducing magnetic current by 1.5, 3
//double Jcu = (Jsc/2.0)/1.5; // current density in copper

double f_he = 0.25; // helium fraction
double dsheath = pow(Btmax,2)/(2.0*muo*sigma_sheath);   // sheath thickness
double f_sheath = 1.0 - pow(1.0-dsheath, 2); // sheath fraction
double f_ins = 0.11; // insulator fraction


// Cross-sectional areas:

Itf = 40000.0;
Acu = Itf/Jcu;
Asc = Itf/Jsc;
Astrands = Acu + Asc;
Aturn = Astrands/(1.0 - f_he - f_sheath - f_ins);
Ains = Aturn*f_ins;
Asheath = Aturn*f_sheath;


// Coil Casing Cross Section

double R2 = 1.03*3.87*rTFinner; // outer external structure radius
int Ntf = 16;  // number of TF coils, make an input
b = 0.02; // side thickness of the casing is b/2
double ATF = Aturn*Bt*2.0*pi*(*Machine).plasma_major_radius()/(Ntf*Itf*muo);
Nwp = int(ATF/Aturn); // number of turns per winding pack


// Initial guess for casing thicknesses

double a, R1;
  a=b;
for(int j=0; j<5; j++)
{


// Calculate R1

R1= rTFinner - 2*a;
for(int i=0; i<5; i++)
{
    R1 = rTFinner - a - ATF/( 2*pi*R1/Ntf - b);
}
wTF = 2*pi*R1/Ntf - b;   // width of the TF coil
dTF = rTFinner - a - R1;  // thickness of the TF coil


// Structure Calculation based on bending:

double G1 = pow(3/(16*sigmab*muo), 0.5);
double G2 = Bt*(*Machine).plasma_major_radius()/(R2-R1);
double G3 = pow( (R2*R2-R1*R1)*(1+log(R2/R1)) + R1*R2*(  2 - (log(R2/R1))*(log(R2/R1))  ),0.5 );
a = 0.208*G1*G2*G3;
}


// Defining parameters of the cross section

delta_inner = a + dTF;
delta_outer = a + dTF;
delta_top = a + dTF;

Btmax = Bt*(*Machine).plasma_major_radius()/(rTFinner - delta_inner/2.0);


// Get the inner-most and outer-most radius of the TF Coil:

  if(kbtmax == 2) (*Machine).get_tf_coil_radii(R1, R2);

} // End of iterative loop over Btmax


double delta = max(delta_inner, max(delta_outer, delta_top));
double width = wTF + b;


// Cross Sectional Area of TF Coil and Volume Fractions

double Acomp = delta*width;    // cross sectional area of the composition
// ... steel
double Asteel = delta*width - dTF*wTF; // cross sectional area of the steel
double vfrac_steel =  Asteel/Acomp;
// ... superconductor
double vfrac_sc = Asc*Nwp/Acomp;
// ... copper
double vfrac_cu = Acu*Nwp/Acomp;
// ... sheath
double vfrac_sheath = Asheath*Nwp/Acomp;
// ... insulator
double vfrac_ins = Ains*Nwp/Acomp;


// Get TF Cross-Section Material Components

   (*Magnet).get_cross_section(0, delta, width);
   (*Magnet).get_material_component("SS", vfrac_steel, 1.0, 0);
   (*Magnet).get_material_component("YBCO", vfrac_sc, 1.0, 1);
   (*Magnet).get_material_component("Cu", vfrac_cu, 1.0, 2);
   (*Magnet).get_material_component("Inconel", vfrac_sheath, 1.0, 3);
   (*Magnet).get_material_component("Polyamide", vfrac_ins, 1.0, 4);
   (*Magnet).get_material_component("Void", 0.0, 1.0, 5);
   
} // eof


// ====================================================================================
// ... r-theta plane, Nb3Sn

void design_tf_coil_rtheta_Nb3Sn(DesignPoint * Machine, 
	Part * Magnet, Part (*RefBody)[2], Part * Plasma, double Bt)
{
if( (*Machine).data_point_filter(0) > 0 ) return;
   const double pi = 4.0*atan(1.0);
   const double muo = 4.0*pi*1.e-7;
   const double sigma_sheath = 8.0e8;
   const double sigmam = 900.0e6;
   const double sigmab = 600.0e6;


// Define the bounding box for the plasma
     double ksiloP = (*Plasma).bounding_box(0, "outer");
     

// Define the bounding box for the Vacuum Vessel
     double ksiloRB = (*RefBody)[0].bounding_box(0, "outer");
     double ksihiRB = (*RefBody)[0].bounding_box(2, "outer");

// Define the Inboard radius of the TF coil that corresponds to its inner
// surface. The inboard thicknesses of all the power core elements enclosed
// within the TF coil including the gaps in between are subtracted from the
// ksiloP of the plasma.

double deltainb = 0.0;
for( int i=0; i<9; i++)
{
     deltainb = deltainb + (*Machine).inboard_generic_build(i);
}
double rTFinner = ksiloP -  deltainb;


// Maximum Magnetic Field

double Btmax = Bt*(*Machine).plasma_major_radius()/rTFinner;


// Current Density in Superconductor versus Btmax function for Nb3Sn @ 4.2K

double btmaxVec[5] = {6.0, 10.0, 14.0, 16.0, 18.0};
double jscVec[5] = {1.52e9, 1.35e9, 8.84e8, 5.43e8, 1.30e8};

// LCC modified jscVec 6/29/10 to determine effect of reducing magnetic current by 1.5, 3
// current estimates may be too optimistic 
//double jscVec[5] = {(1.52e9)/1.5, (1.35e9)/1.5, (8.84e8)/1.5, (5.43e8)/1.5, (1.30e8)/1.5};

double Itf, Aturn, Astrands, Acu, Asc, Asheath, Ains;
double delta_inner, delta_outer, delta_top;
double wTF, dTF, b;
int Nwp;

for(int kbtmax=0; kbtmax<3; kbtmax++)  // Iterative loop over Btmax
{
if(Btmax < btmaxVec[0] || Btmax > btmaxVec[4])
{
  cout << "Btmax out of range\n";         // LCC, this error is seen sometimes in sysoutfinal log
  (*Machine).get_data_point_filter(1, 1);
  return;
}
double dist = fabs(Btmax - btmaxVec[0]);
int imark = 0;

for(int i=0; i<5; i++)
{
if( fabs(Btmax - btmaxVec[i]) < dist)
{
  dist = fabs(Btmax - btmaxVec[i]);
  imark = i;
}
}


// Cross section of the winding pack:

double Jsc = jscVec[imark]; // current density of the superconductor
double Jcu = 2.24e8; // current density in copper

// LCC modified Jcu 6/29/10 to determine effect of reducing magnetic current by 1.5, 3
//double Jcu = (2.24e8)/1.5; // current density in copper

double f_he = 0.25; // helium fraction
double dsheath = pow(Btmax,2)/(2.0*muo*sigma_sheath);   // sheath thickness
double f_sheath = 1.0 - pow(1.0-dsheath, 2); // sheath fraction
double f_ins = 0.1; // insulation


// Cross-sectional areas:

Itf = 40000.0;
Acu = Itf/Jcu;
Asc = Itf/Jsc;
Astrands = Acu + Asc;
Aturn = Astrands/(1.0 - f_he - f_sheath - f_ins);
Ains = Aturn*f_ins;
Asheath = Aturn*f_sheath;


// Coil Casing Cross Section

double A = pow( Bt*(*Machine).plasma_major_radius(), 2)/(2.0*muo);
double R2 = rTFinner + (ksihiRB-ksiloRB)*1.3; // outer external structure radius
int Ntf = 16;  // number of TF coils
b = 0.02; // side thickness of the casing is b/2
double ATF = Aturn*Bt*2.0*pi*(*Machine).plasma_major_radius()/(Ntf*Itf*muo);
Nwp = int(ATF/Aturn); // number of turns per winding pack


// Initial guess for casing thicknesses

double t1 = b;
double t2 = b;
double a = b;
double R1, R3;
double S1, S2;

for(int j=0; j<5; j++)
{

// Calculate R1

R1= rTFinner - 2*t1;
for(int i=0; i<5; i++)
{
    R1 = rTFinner - t1 - ATF/( 2*pi*R1/Ntf - b);
}

wTF = 2*pi*R1/Ntf - b;     // width of the TF coil
dTF = rTFinner - t1 - R1;  // thickness of the TF coil


// Structure Calculation based on bending:

S1 = A*( R2*log(R2/R1) - R2 + R1 )/( R1*(R2-R1)  );
S2 = A*( R2 - R1 - R1*log(R2/R1)  )/( R2*(R2-R1)  );
t1 = S1/sigmam;
t2 = 1.5*S2/sigmam;
R3 = ( S1*R1*R1 + S2*R2*R2  )/(A*log(R2/R1));

// a = pow( 3.0/8.0*S2*R2*(R2-R3)/(sigmab*R3), 0.5  );
// a with "correction"
// previous trial: a = pow( 6.0*S2*R2*(R2-R3)/(sigmab*R3), 0.5  );

double G1 = pow(3.0/(16.0*sigmab*muo), 0.5);
double G2 = Bt*(*Machine).plasma_major_radius()/(R2-R1);
double G3 = pow( (R2*R2-R1*R1)*(1+log(R2/R1)) + R1*R2*(  2 - (log(R2/R1))*(log(R2/R1))  ),0.5 );
  a = 0.3551*G1*G2*G3;
}


// Defining parameters of the cross section

delta_inner = t1 + dTF;
delta_outer = t2 + dTF;
delta_top = a + dTF;

Btmax = Bt*(*Machine).plasma_major_radius()/(rTFinner - delta_inner/2.0);


} // End of iterative loop over Btmax

double delta = max(delta_inner, max(delta_outer, delta_top));
double width = wTF + b;


// Cross Sectional Area of TF Coil and Volume Fractions

double Acomp = delta*width;    // cross sectional area of the composition
// ... steel
double Asteel = delta*width - dTF*wTF; // cross sectional area of the steel
double vfrac_steel =  Asteel/Acomp;
// ... superconductor
double vfrac_sc = Asc*Nwp/Acomp;
// ... copper
double vfrac_cu = Acu*Nwp/Acomp;
// ... sheath
double vfrac_sheath = Asheath*Nwp/Acomp;
// ... insulator
double vfrac_ins = Ains*Nwp/Acomp;


// Get TF Cross-Section Material Components

   (*Magnet).get_cross_section(0, delta, width);
   (*Magnet).get_material_component("SS", vfrac_steel, 1.0, 0);
   (*Magnet).get_material_component("Nb3Sn", vfrac_sc, 1.0, 1);
   (*Magnet).get_material_component("Cu", vfrac_cu, 1.0, 2);
   (*Magnet).get_material_component("Inconel", vfrac_sheath, 1.0, 3);
   (*Magnet).get_material_component("Polyamide", vfrac_ins, 1.0, 4);
   (*Magnet).get_material_component("Void", 0.0, 1.0, 5);
   
} // eof


// ====================================================================================
// ... r-theta plane, ARIES-AT coil design with fixed volume fractions ...

void design_tf_coil_rtheta_aries_at(DesignPoint * Machine, Part * Magnet, Part * Plasma, double Bt)
{
if( (*Machine).data_point_filter(0) > 0 ) return;
   const double pi 		= 4.0*atan(1.0);
   const double muo 	= 4.0*pi*1.e-7;
   const double sigmam 	= 900.0e6;
   const double sigmab 	= 600.0e6;
   

// Define the bounding box for the plasma
     double ksiloP = (*Plasma).bounding_box(0, "outer");

// Define the Inboard radius of the TF coil that corresponds to its inner
// surface. The inboard thicknesses of all the power core elements enclosed
// within the TF coil including the gaps in between are subtracted from the
// ksiloP of the plasma.

double deltainb = 0.0;
for( int i=0; i<9; i++)
{
     deltainb = deltainb + (*Machine).inboard_generic_build(i);
}
double rTFinner = ksiloP -  deltainb;


// Maximum Magnetic Field
double Btmax = Bt*(*Machine).plasma_major_radius()/rTFinner;


// Current Density in Superconductor versus Btmax function for YBCO @ 75K
double btmaxVec[5] = {6.0, 10.0, 14.0, 16.0, 18.0};
double jscVec[5] =   {9.00e8, 6.75e8, 5.40e8, 4.86e8, 4.62e8};
double Itf, Aturn, Aag, Asc, Asheath, Aceo2, Apolyamide;
double delta_inner, delta_outer, delta_top;
double wTF, dTF, b;
int Nwp;

for(int kbtmax=0; kbtmax<3; kbtmax++)  // Iterative loop over Btmax
{
if(Btmax < btmaxVec[0] || Btmax > btmaxVec[4])
{
  cout << "Btmax out of range\n";
  (*Machine).get_data_point_filter(1, 1);
  return;
}
double dist = fabs(Btmax - btmaxVec[0]);
int imark = 0;

for(int i=0; i<5; i++)
{
if( fabs(Btmax - btmaxVec[i]) < dist)
{
  dist = fabs(Btmax - btmaxVec[i]);
  imark = i;
}
}


// Cross section of the winding pack:
double Jsc = jscVec[imark]; // current density of the superconductor


// Cross-sectional areas:
// ... local volume fractions ...
//    72% Inconel-625, 0.5% Ag, 7% YBa2Cu3O5, 7% CeO2 insulator, 13.5% GFF polyimide insulator

Itf = 40000.0;
double f_sc = 0.07;
double f_ag = 0.005;
double f_sheath = 0.72;
double f_CeO2 = 0.07;
double f_polyamide = 0.135;

Asc = Itf/Jsc;
Aturn =  Asc/f_sc;
Aag = Aturn*f_ag;
Asheath = Aturn*f_sheath;
Aceo2 = Aturn*f_CeO2;
Apolyamide = Aturn*f_polyamide;


// Coil Casing Cross Section

double A = pow( Bt*(*Machine).plasma_major_radius(), 2)/(2.0*muo);
double R2 = 3.87*rTFinner; // outer external structure radius
int Ntf = 16;  // number of TF coils
b = 0.02; // side thickness of the casing is b/2
double ATF = Aturn*Bt*2.0*pi*(*Machine).plasma_major_radius()/(Ntf*Itf*muo);
Nwp = int(ATF/Aturn); // number of turns per winding pack


// Initial guess for casing thicknesses

double t1 = b;
double t2 = b;
double a = b;
double R1, R3;
double S1, S2;

for(int j=0; j<5; j++)
{


// Calculate R1

R1= rTFinner - 2*t1;
for(int i=0; i<5; i++)
{
    R1 = rTFinner - t1 - ATF/( 2*pi*R1/Ntf - b);
}

wTF = 2*pi*R1/Ntf - b;   // width of the TF coil
dTF = rTFinner - t1 - R1;  // thickness of the TF coil


// Structure Calculation based on bending:

S1 = A*( R2*log(R2/R1) - R2 + R1 )/( R1*(R2-R1)  );
S2 = A*( R2 - R1 - R1*log(R2/R1)  )/( R2*(R2-R1)  );
t1 = S1/sigmam;
t2 = 1.5*S2/sigmam;
R3 = ( S1*R1*R1 + S2*R2*R2  )/(A*log(R2/R1));
a = pow( 3.0/8.0*S2*R2*(R2-R3)/(sigmab*R3), 0.5  );
}


// Defining parameters of the cross section

delta_inner = t1 + dTF;
delta_outer = t2 + dTF;
delta_top = a + dTF;

Btmax = Bt*(*Machine).plasma_major_radius()/(rTFinner - delta_inner/2.0);

} // End of iterative loop over Btmax


double delta = max(delta_inner, max(delta_outer, delta_top));
double width = wTF + b;


// Cross Sectional Area of TF Coil and Volume Fractions

double Acomp = delta*width;    // cross sectional area of the composition
// ... steel
double Asteel = (delta*width - dTF*wTF)*0.95; // cross sectional area of the steel
// double Aln = Asteel*0.05/0.95; // cross sectional area of the liquid nitrogen
// ... steel
double vfrac_steel =  Asteel/Acomp;
// ... superconductor
double vfrac_sc = Asc*Nwp/Acomp;
// ... copper
double vfrac_ag = Aag*Nwp/Acomp;
// ... sheath
double vfrac_sheath = Asheath*Nwp/Acomp;
// ... CeO2 insulation ....
double vfrac_ceo2 = Aceo2*Nwp/Acomp;
// ... polyamide insulation
double vfrac_polyamide = Apolyamide*Nwp/Acomp;


// Get TF Cross-Section Material Components

   (*Magnet).get_cross_section(0, delta, width);
   (*Magnet).get_material_component("SS", vfrac_steel, 1.0, 0);
   (*Magnet).get_material_component("YBCO", vfrac_sc, 1.0, 1);
   (*Magnet).get_material_component("Ag", vfrac_ag, 1.0, 2);
   (*Magnet).get_material_component("Inconel", vfrac_sheath, 1.0, 3);
   (*Magnet).get_material_component("CeO2", vfrac_ceo2, 1.0, 4);
   (*Magnet).get_material_component("Polyamide", vfrac_polyamide, 1.0, 5);

} // eof


// ====================================================================================
// ********* CURRENT TF r-z algorithm in use ******************************************
// ====================================================================================
// ... r-z plane (revolved)

void design_tf_coil_rz_plane(DesignPoint * Machine, Part * Magnet, Part (*VVRefBody)[4])
 									// VVRefBody = VacuumVesselDCLL (DCLL)
 									// VVRefBody = VacuumVessel     (SiC)
{													
if( (*Machine).data_point_filter(0) > 0 ) return;

// Vacuum Vessel Box:
double ksihi, etahi;
if (blanketmode == "dcll") {
       ksihi = (*VVRefBody)[1].bounding_box(2, "outer");   // for DCLL, [1] is OB
       etahi = (*VVRefBody)[2].bounding_box(3, "outer"); } // for DCLL, [2] is top
else { ksihi = (*VVRefBody)[0].bounding_box(2, "outer");   // for SiC,  [0] is "main" VV
       etahi = (*VVRefBody)[0].bounding_box(3, "outer"); }

// Wall and wall gap thickness
    double dw    = (*Magnet).cross_section_width(0);   // TF magnet wall thickness
//  double dwgap = (*Machine).inboard_generic_build(17); // VV and TF magnet gap, not in use

// X-point height
    double hxp   = (*Machine).xpoint(1);
    
// TF Coil inboard and outboard radii:
    double R1 = (*Machine).tf_magnet_inner_leg_radius();
    double R2 = (*Machine).tf_magnet_outer_leg_radius();
    double r1 = R1 + dw;
    double r2 = R2 - dw;
    
// Plasma major radius:
   double Ro = (*Machine).plasma_major_radius();
   double ao = (*Machine).plasma_minor_radius();
    
// Define the parameters of the "Princeton-D" ellipse:
// ...from Richard J. Thome, John M. Tarrh, ÒMHD and Fusion Magnets Ð Field and Force Design ConceptsÓ, John Wiley & Sons, 1982.
   double b  = etahi + 1.043 ; // radius on the vertical axis, same as cap below, reference?
   double a1 = (Ro-r1)/pow(1.0-pow( 0.95*etahi/b, 2.0 ), 0.5 );
   double a2 = r2 - Ro;

// Transform the plasma into the TF coil contour
// LCC - does this mean the plasma contour is stretched to the D-ellipse shape???
  int len = 300;
  double *x, *y;
  x = new double [len];
  y = new double [len];
  
  x[0] = r1;
  y[0] = 0.0;
  double delx = 2*(r2-r1)/(len-3);
  for(int i=1; i<len-1; i++)		// loop over TF coil contour
  { if(i < (len+1)/2 ) 				// upper half of the coil
	{ int c;
      if(i == 1){c=0;}
      else      {c=1;}
      	x[i]=x[i-1]+c*delx;
      if(x[i] < (r2-a2)) 			// inner ellipse
        { y[i] = b*pow(1.0-pow((x[i]-Ro)/a1, 2.0), 0.5); }
      else  						// outer ellipse
        { y[i] = b*pow(1.0-pow((x[i]-Ro)/a2, 2.0), 0.5); } }
      else   						// lower half of the coil
        { if(i == (len+1)/2) 
        { x[i] = r2; }
      else
        { x[i] = x[len-i]; }
      if(i == len-2) x[i] = x[0];
      if(x[i] < (r2-a2)) 			// inner ellipse
        { y[i] = -b*pow(1.0-pow((x[i]-Ro)/a1, 2.0), 0.5); }
      else    						// outer ellipse
        { y[i] = -b*pow(1.0-pow((x[i]-Ro)/a2, 2.0), 0.5); } }
  } // end loop over TF coil contour
  
// Final connection point for inner loop
      x[len-1] = r1;
      y[len-1] = -delx;


// Update the Magnet
// ... inner loop ...
	(*Magnet).get_contour_size(len, "inner");
	(*Magnet).get_contour_size(len, "outer");

	double (*c_in)[2], (*c_out)[2], * d_in, * d_out;

	c_in  = new double [len][2];
	c_out = new double [len][2];
	d_in  = new double [len];
	d_out = new double [len];

	for(int i=0; i<len; i++)
	{ c_in[i][0] = x[i];
 	  c_in[i][1] = y[i];
 	  d_in[i] = (*Magnet).cross_section_depth(0); }


// Define the outer contour at prescribed distance from the inner contour
	c_out[0][0] = c_in[0][0] - dw; // dw = thickness of TF coil
    c_out[0][1] = 0.0;

// Construct the outer curve by looping along the inner curve
for(int i=1; i<len; i++) 
	{ double tagm = sqrt( pow( (c_in[i+1][0]-c_in[i-1][0]),2) +
                          pow( (c_in[i+1][1]-c_in[i-1][1]),2) );
	  double nx = - (c_in[i+1][1]-c_in[i-1][1])/tagm;
	  double ny =   (c_in[i+1][0]-c_in[i-1][0])/tagm;
	
	if( nx < -0.9 ) // slopes of ellipses
	{ c_out[i][0] = c_in[i][0] + nx*dw;
	  c_out[i][1] = c_in[i][1] + ny*dw; }

	if(nx >= -0.9 && nx < -0.7)
	{ c_out[i][0] = c_in[i][0] + nx*( dw*(nx+0.9)/0.2 + dw*(nx+0.7)/(-0.2) );
	  c_out[i][1] = c_in[i][1] + ny*( dw*(nx+0.9)/0.2 + dw*(nx+0.7)/(-0.2) ); }
    
    if(nx >= -0.7)
	{ c_out[i][0] = c_in[i][0] + nx*dw;
	  c_out[i][1] = c_in[i][1] + ny*dw;}
} // end loop along the inner curve 

// Final connection point for outer loop
     c_out[len-1][0] = c_in[0][0] - dw;
     c_out[len-1][1] = 0.0;


for(int i=0; i<len; i++)
{ d_out[i] =   (*Magnet).cross_section_depth(0); } // 1.03758
(*Magnet).get_contours(c_in, c_out, d_in, d_out);


// Volume and Surface Areas
(*Magnet).get_volume();
(*Magnet).get_surface_areas();
(*Magnet).get_perimeters();

delete [] x;
delete [] y;
delete [] c_in;
delete [] c_out;
delete [] d_in;
delete [] d_out;

} // eof


// *********** End of coil designs ****************************************************

// ====================================================================================
// Bucking Cylinder

void design_bucking_cylinder(DesignPoint * Machine, Part * BuckCyl, 
							 Part * Magnet, 		Part (*Vessel)[4], double Bt)
{ 
if( (*Machine).data_point_filter(0) > 0 ) return;
   const double pi  = 4.0*atan(1.0);
   const double muo = 4.0*pi*1.e-7;

// Open the input file to read the magnet parameters:
  ifstream fin_mag;  fin_mag.open("InputData/magnets.data");
  fin_mag.width(10); fin_mag.precision(2);
  fin_mag.setf(ios_base::uppercase, ios_base::scientific);
  int itxt_mag = 57, inum_mag = 13;
  int iskip_mag = itxt_mag + inum_mag;

// ... Yield stress for the bucking cylinder
   double sigmaYield; fin_mag.seekg(8*iskip_mag + itxt_mag + 1); fin_mag >> sigmaYield;
   double sigma = sigmaYield;

// ... Define the Inner Contour ...
// ... start from the outer bounding box of TF coil as a reference entity and build on top of it

// Define the bounding box of the Vacuum Vessel, top of VV
	double etahiV;
if (blanketmode == "dcll") 
     etahiV = (*Vessel)[2].bounding_box(3, "outer");
else etahiV = (*Vessel)[0].bounding_box(3, "outer");

// ... ... reference points ... ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];

     double ksilo =(*Magnet).bounding_box(0, "outer");
//     double ksilo =(*Magnet).bounding_box(0, "inner");   // altered 7/10/12 by M. Tillack
     double dscr  = (*Machine).inboard_generic_build(21);  // TF coil to bucking cyl gap

     double rbtmax =(*Magnet).bounding_box(0, "inner");
     double Btmax = Bt*(*Machine).plasma_major_radius()/rbtmax;

// original calculation of bucking cylinder thickness
//   double dwall = pow(Btmax,2)*ksilo/(2.0*muo*sigma);

// modified calculation of dwall (7/10/12) by M. Tillack
       double rtfin  = (*Magnet).bounding_box(0, "outer");
       double rtfout = (*Magnet).bounding_box(0, "inner");
       double rtfave = (rtfin + rtfout)/2.0;
       double rbcout = rtfin - dscr;
       double press = (rbcout/rtfave) * (pow(Btmax,2)/(2.*muo));
       double rbcin = rbcout * sqrt(1.-(2.*press)/(sigma));
       double dwall = rbcout - rbcin;
//   rtfout is the outer TF coil radius on IB
//    rtfin is the inner TF coil radius on IB (rtfin < rtfout)
//   rbcout is the outer bucking cylinder radius on IB
//    rbcin is the inner bucking cylinder radius (rbcin < rbcout)
//    dwall is the bucking cylinder thickness

// Define the Outer Contour
// ... reference points ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 4;
     double xs[10], ys[10];

     xs[0] = ksilo - dscr - dwall;
     ys[0] = 0.0;

     xs[1] = xs[0];
     ys[1] = 0.95*etahiV;

     xs[2] = ksilo - dscr;
     ys[2] = ys[1];

     xs[3] = xs[2];
     ys[3] = 0.0;

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

// ... Close the loop by adding a symmetric profile in the clockwise direction: ...
     for(int i=0;  i<nptso; i++)
     { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
       ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
       nptso = 2*nptso-2;

// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line"; }

// Update the Bucking Cylinder
(*BuckCyl).get_contour_descriptions(nptsi,        nptso,		resol_i,      resol_o,
                                    ref_pts_i,    ref_pts_o,	connectors_i, connectors_o);


// Define the Inner Target Curve For the Bucking Cylinder
double targcurve[10000][2]; int tgl=0;


// Connect the Reference Points into Contours
(*BuckCyl).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation
 double * d_in; double * d_out;
 int lin  = 1;
 int lout = (*BuckCyl).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*BuckCyl).contour(i, 0, "outer")/16.0;  // LCC part is revolved here, /16 for unit
 (*BuckCyl).get_depth(d_in, d_out);


// Volume and Surface Areas calculated in Part.cpp
(*BuckCyl).get_volume();
(*BuckCyl).get_surface_areas();

delete [] d_in;
delete [] d_out;

fin_mag.close();
} // eof



// ====================================================================================
// Upper Toroidal Cap on the TF Coil

void design_upper_cap_shape(DesignPoint * Machine, 	Part * UpperCap,
							Part * RefBody, 		Part (*Vessel)[4], double Bt)
							// RefBody = TF coil
 { if( (*Machine).data_point_filter(0) > 0 ) return;

// Constants:
  const double pi    = 4.0*atan(1.0); // pi = 3.14159....
  const double muo   = 4.0*pi*1.e-7;

// Open the input file to read the magnet parameters:
  ifstream fin_mag;  fin_mag.open("InputData/magnets.data");
  fin_mag.width(10); fin_mag.precision(2);
  fin_mag.setf(ios_base::uppercase, ios_base::scientific);
  int itxt_mag = 57, inum_mag = 13;
  int iskip_mag = itxt_mag + inum_mag;

// ... Yield stress for the bucking cylinder
   double sigmaYield; fin_mag.seekg(8*iskip_mag + itxt_mag + 1); fin_mag >> sigmaYield;
   double sigma = 2.0/3.0*sigmaYield;

// ... Define the Inner Contour ...
     double ksilo =(*RefBody).bounding_box(0, "outer");
//   double ksihi =(*RefBody).bounding_box(2, "outer"); // LCC turn this on to 2~hi?
//   double etahi =(*RefBody).bounding_box(3, "outer"); 
     

// Define the bounding box of the Vacuum Vessel, top of VV
	double etahiV;
if (blanketmode == "dcll")
	 etahiV = (*Vessel)[2].bounding_box(3, "outer");
else etahiV = (*Vessel)[0].bounding_box(3, "outer");
//cout << "etahiV = " << (*Vessel)[0].bounding_box(3, "outer") << "\n";     

// Toroidal Cap thickness and gap from the TF coil
     double Btmax   = Bt*(*Machine).plasma_major_radius()/ksilo;
     double dwgap   = (*Machine).vertical_generic_build(18);  // TF coil to toroidal cap gap
     double dwall   = pow(Btmax,2)*ksilo/(2.0*muo*sigma);     // toroidal cap thickness
     double dmagnet = (*RefBody).cross_section_width(0);      // TF magnet thickness
//cout << "dwall 1 = " << dwall << "\n";     

// ... ... reference points ... ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];


// Define the Outer Contour
// ... ... reference points ... ...
   double ref_pts_o[100][2], center_o[100][2];
   int nptso = 4;
   double Ro = (*Machine).plasma_major_radius();
   double R1 = (*Machine).tf_magnet_inner_leg_radius();
   double R2 = (*Machine).tf_magnet_outer_leg_radius();
   double r1 = R1 + dmagnet;
   double r2 = R2 - dmagnet;
   double b = etahiV + 1.043; // TF ellipse radius on the vertical axis, -etaloV in lower
   double a2 = r2 - Ro;       // TF outer ellipse radius on the horizontal axis


// Limiting point on the outer TF coil contour:
// ... this point is determined by the height of the vacuum vessel etahiV ...
    double y_tf = etahiV;
    double x_tf = Ro + a2*sqrt(1.0 - pow(y_tf/b,2.0));
    double g = (x_tf - Ro)/a2;
    double alpha_t = atan( -b/a2*g/sqrt(1.0 - pow(g,2.0)) ); // LCC +b in lowercap
    double alpha_n = alpha_t + pi/2.0; // LCC -pi in lowercap
    double xs[10], ys[10];
    
// LCC different order [0]-[3] in lowercap
      xs[0] = ksilo - dwall - dwgap;
      ys[0] = 0.95 * etahiV + dwgap; // 95% of height of VV top

      xs[1] = x_tf + (dwgap + dwall + dmagnet)*cos(alpha_n);
      ys[1] = y_tf + (dwgap + dwall + dmagnet)*sin(alpha_n);

      xs[2] = x_tf + (dwgap + dmagnet)*cos(alpha_n);
      ys[2] = y_tf + (dwgap + dmagnet)*sin(alpha_n);

      xs[3] = xs[0] + dwall;
      ys[3] = ys[0];
    
/*cout << "xs[0], ys[0] = " << xs[0] << ",	" << ys[0] << "\n";
cout << "xs[1], ys[1] = " << xs[1] << ",	" << ys[1] << "\n";
cout << "xs[2], ys[2] = " << xs[2] << ",	" << ys[2] << "\n";
cout << "xs[3], ys[3] = " << xs[3] << ",	" << ys[3] << "\n";
*/
     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
      if(i == 0 || i == 2) connectors_o[i] = "target"; }


// Update the Blanket
(*UpperCap).get_contour_descriptions(nptsi, 	nptso, 		resol_i, 		resol_o,
                                 	 ref_pts_i, ref_pts_o, 	connectors_i, 	connectors_o);


// Define the Inner Target Curve For Blanket Shape
   double targcurve[10000][2];
   int tgl = (*RefBody).contour_size("outer");
    for(int i=0; i<tgl; i++)
    { targcurve[i][0] = (*RefBody).contour(i, 0, "outer");
      targcurve[i][1] = (*RefBody).contour(i, 1, "outer"); }


// Connect the Reference Points into Contours
(*UpperCap).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*UpperCap).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*UpperCap).contour(i, 0, "outer")/16.0;
 	(*UpperCap).get_depth(d_in, d_out);

// Volume and Surface Areas
	(*UpperCap).get_volume();
	(*UpperCap).get_surface_areas();
//cout << "upper cap vol = " << (*UpperCap).volume()*16.0 << "\n";

delete [] d_in;
delete [] d_out;
} // eof


// ====================================================================================
// Lower Toroidal Cap on the TF Coil - By Design (not in use now)

// LCC NEED to fix vol of lower cap, should be = to upper cap. Try following ReplHT shield algorithm?
// DCLL module uses  copy_upper_cap_into_lower_cap_shape_by_symmetry(&LowerCap, &UpperCap);

void design_lower_cap_shape(DesignPoint * Machine, 	Part * LowerCap,
							Part * RefBody, 		Part (*Vessel)[4], double Bt)

 { if( (*Machine).data_point_filter(0) > 0 ) return;

// Constants:
  const double pi    = 4.0*atan(1.0); // pi = 3.14159....
  const double muo   = 4.0*pi*1.e-7;

// Open the input file to read the magnet parameters:
  ifstream fin_mag;  fin_mag.open("InputData/magnets.data");
  fin_mag.width(10); fin_mag.precision(2);
  fin_mag.setf(ios_base::uppercase, ios_base::scientific);
  int itxt_mag = 57, inum_mag = 13;
  int iskip_mag = itxt_mag + inum_mag;

// ... Yield stress for the bucking cylinder
   double sigmaYield; fin_mag.seekg(8*iskip_mag + itxt_mag + 1); fin_mag >> sigmaYield;
   double sigma = 2.0/3.0*sigmaYield;
   
// ... Define the Inner Contour ...
/* ... ... start from the first wall as the reference entity and build on top of it ... */

     double ksilo =(*RefBody).bounding_box(0, "outer");
//     double ksihi =(*RefBody).bounding_box(2, "outer");
//     double etahi =(*RefBody).bounding_box(3, "outer");


// Define the bounding box of the Vacuum Vessel
	double etaloV;
if (blanketmode == "dcll")
	 etaloV = (*Vessel)[3].bounding_box(1, "outer"); // maybe 3, "outer"?
else etaloV = (*Vessel)[0].bounding_box(3, "outer");


// Toroidal Cap thickness and gap from the TF coil
     double Btmax   = Bt*(*Machine).plasma_major_radius()/ksilo;
     double dwgap   = (*Machine).vertical_generic_build(18); // TF coil to upper cap gap
     double dwall   = pow(Btmax,2)*ksilo/(2.0*muo*sigma);    // toroidal cap thickness
     double dmagnet = (*RefBody).cross_section_width(0);     // wall thickness of TF magnet
//cout << "dwgap = " << dwgap  << "\n";     

// ... ... reference points ... ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];

// Define the Outer Contour
// ... ... reference points ... ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 4;
     double Ro = (*Machine).plasma_major_radius();
     double R1 = (*Machine).tf_magnet_inner_leg_radius();
	 double R2 = (*Machine).tf_magnet_outer_leg_radius();
	 double r1 = R1 + dmagnet;
	 double r2 = R2 - dmagnet;
	 double b = -etaloV + 1.043; // TF ellipse radius on the vertical axis
	 double a2 = r2 - Ro;        // TF outer ellipse radius on the horizontal axis


// Limiting point on the outer TF coil contour:
// ... this point is determined by the height of the vacuum vessel etahiV ...
    double y_tf = etaloV;
    double x_tf = Ro + a2*sqrt(1.0 - pow(y_tf/b,2.0));
    double g = (x_tf - Ro)/a2;
    double alpha_t = atan( b/a2*g/sqrt(1.0 - pow(g,2.0)) ); // -b in uppercap
    double alpha_n = alpha_t - pi/2.0; // +pi in uppercap
    double xs[10], ys[10];
      
      xs[0] = ksilo - dwall - dwgap;
      ys[0] = 0.95*etaloV - dwgap;
      
      xs[1] = xs[0] + dwall;
      ys[1] = ys[0];
      
      xs[2] = x_tf + (dwgap + dmagnet)*cos(alpha_n);
      ys[2] = y_tf + (dwgap + dmagnet)*sin(alpha_n);

      xs[3] = x_tf + (dwgap + dwall + dmagnet)*cos(alpha_n);
      ys[3] = y_tf + (dwgap + dwall + dmagnet)*sin(alpha_n);

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }

// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
    if(i == 1 || i == 3) connectors_o[i] = "target"; }


// Update the Blanket
(*LowerCap).get_contour_descriptions(nptsi, nptso, resol_i, resol_o, ref_pts_i, 
									 ref_pts_o, connectors_i, connectors_o);

// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2];
int tgl = (*RefBody).contour_size("outer");
    for(int i=0; i<tgl; i++)
    { targcurve[i][0] = (*RefBody).contour(i, 0, "outer");
      targcurve[i][1] = (*RefBody).contour(i, 1, "outer"); }


// Connect the Reference Points into Contours
(*LowerCap).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*LowerCap).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*LowerCap).contour(i, 0, "outer")/16.0;
 	(*LowerCap).get_depth(d_in, d_out);

// Volume and Surface Areas
	(*LowerCap).get_volume();
	(*LowerCap).get_surface_areas();
//cout << "lower cap vol = " << (*LowerCap).volume()*16.0 << "\n"; // should be ~73 like upper

delete [] d_in;
delete [] d_out;
} // eof


// ====================================================================================
// Lower Toroidal Cap on the TF Coil - By Symmetry (from Upper Cap)
// LCC - USE this rather than "by design" above

void copy_upper_cap_into_lower_cap_shape_by_symmetry(Part * LowerCap, Part* UpperCap)
{

 int lin  = (*UpperCap).contour_size("inner");
 int lout = (*UpperCap).contour_size("outer");

 const double pi = 4.0*atan(1.0);
 double * d_in; 
 double * d_out;
  d_in  = new double [lin];
  d_out = new double [lout];

(*LowerCap).get_contour_size(lin, "inner");
(*LowerCap).get_contour_size(lout, "outer");

double (*ctr_in)[2], (*ctr_out)[2];

ctr_in  = new double [lin][2];
ctr_out = new double [lout][2];

for(int i=0; i<lin; i++)
{ ctr_in[i][0] = 0.0;
  ctr_in[i][1] = 0.0;
  d_in[i] = 0.0; }

for(int i=0; i<lout; i++)
{ ctr_out[i][0] =  (*UpperCap).contour(lout-1-i, 0, "outer");
  ctr_out[i][1] = -(*UpperCap).contour(lout-1-i, 1, "outer");
  d_out[i] = 2.0*pi*ctr_out[i][0]/16.0; }

(*LowerCap).get_contours(ctr_in, ctr_out, d_in, d_out);


// ... Volume and Surface Areas ...
(*LowerCap).get_volume();
(*LowerCap).get_surface_areas();

delete [] d_in;
delete [] d_out;
delete [] ctr_in;
delete [] ctr_out;

} // eof


// ====================================================================================
// Maintenance Port Window inbetween the TF Coils

void design_mpw(DesignPoint * Machine, Part * Window, Part * Magnet, Part(*Vessel)[4])
{
if( (*Machine).data_point_filter(0) > 0 ) return;
     const double pi = 4.0*atan(1.0);


// Define the bounding box for the vacuum vessel, should be same for DCLL as SiC
double etalo, ksihi, etahi;
	etalo = (*Vessel)[0].bounding_box(1, "outer");
	ksihi = (*Vessel)[0].bounding_box(2, "outer");
	etahi = (*Vessel)[0].bounding_box(3, "outer");

     int len = (*Magnet).contour_size("inner");
     double *x, *y;
     x = new double [len+100];
     y = new double [len+100];
     
     
// Define the INNER contour of the maintenance port window
// dimensions???  width must be dis between TF coils at radius of OB VV
     int lwi = 0;
     for(int i=len-1; i>=0; i--)
     { double xmag = (*Magnet).contour(i, 0, "inner");
       double ymag = (*Magnet).contour(i, 1, "inner");
    
     if(xmag > ksihi && ymag >= etalo && ymag <= etahi)
    	{ x[lwi] = xmag;
          y[lwi] = ymag;
          lwi++; } }


// Define the OUTER contour of the maintenance port window
     int lwo = lwi + 50;
     for(int i=0; i<len; i++)
     { double xmag = (*Magnet).contour(i, 0, "outer");
       double ymag = (*Magnet).contour(i, 1, "outer");

     if(xmag > ksihi && ymag >= etalo && ymag <= etahi)
     	{ x[lwo] = xmag;
          y[lwo] = ymag;
          lwo++; } }


// Top contour of the maintenance port window
  for(int i=lwi+1; i<lwi+50; i++)
  	{ x[i] = x[lwi] + (i-lwi)*(x[lwi+50]-x[lwi])/50;
      y[i] = y[lwi]; }


// Bottom contour of the maintenance port window
  for(int i=lwo+1; i<lwo+50; i++)
  	{ x[i] = x[lwo] - (i-lwo)*(x[0]-x[lwo])/50;
      y[i] = y[lwi]; }


// Update the maintenance port window
	int lw = lwo + 50;
	(*Window).get_contour_size(lw, "inner");
	(*Window).get_contour_size(1, "outer");
	double (*c_in)[2], (*c_out)[2], * d_in, * d_out;
	c_in  = new double [lw][2];
	c_out = new double [1][2];
	d_in  = new double [lw];
	d_out = new double [1];

for(int i=0; i<lw; i++)
	{ c_in[i][0]  = x[i];
      c_in[i][1]  = y[i];
  	  d_in[i]     = 2.0*pi*(c_in[i][0]-0.5*ksihi)/16.0; }
  	  c_out[0][0] = 0.0;
  	  c_out[0][1] = 0.0;
  	  d_out[0]    = 0.0;

(*Window).get_contours(c_in, c_out, d_in, d_out);
(*Window).get_volume();
(*Window).get_surface_areas();
  
delete [] x;
delete [] y;
delete [] c_in;
delete [] c_out;
delete [] d_in;
delete [] d_out;

} // eof



// ====================================================================================
// ====================================================================================
// PF Coil Design
// 0-3 top + 32-35 bottom are central solenoid, 4-17 top, 18-31 bottom (CW), 4-17 are spares (top).

void design_pf_coil(DesignPoint * Machine, Part (*Magnet)[36], Part * RefCyl,
						   Part * TopTCap, Part * BottomTCap) // "RefCyl" = Bucking Cylinder
{ if( (*Machine).data_point_filter(0) > 0) return;
    const double pi = 4.0*atan(1.0);


// Open the input file to read the magnet parameters:
  ifstream fin_mag;  fin_mag.open("InputData/magnets.data");
  fin_mag.width(10); fin_mag.precision(2);
  fin_mag.setf(ios_base::uppercase, ios_base::scientific);
  int itxt_mag = 57, inum_mag = 13;
  int iskip_mag = itxt_mag + inum_mag;


// ... Yield stress for the PF coil case [Pa]
   double sigmaYield; fin_mag.seekg(8*iskip_mag + itxt_mag + 1); fin_mag >> sigmaYield;
   const double sigma = 2.0/3.0*sigmaYield; // take 2/3 of yield stress

     
// Define the bounding box of the Bucking Cylinder
     double ksiloBC = (*RefCyl).bounding_box(0, "outer");
//   double etaloBC = (*RefCyl).bounding_box(1, "outer");
//   double ksihiBC = (*RefCyl).bounding_box(2, "outer");
     double etahiBC = (*RefCyl).bounding_box(3, "outer");
//cout << "ksiloBC = " << ksiloBC << "\n";
//cout << "ksiloBC = " << etahiBC << "\n";

// Define the bounding box of the Top Toroidal Cap
     double ksiloT = (*TopTCap).bounding_box(0, "outer");
//   double etaloT = (*TopTCap).bounding_box(1, "outer");
     double ksihiT = (*TopTCap).bounding_box(2, "outer");
//   double etahiT = (*TopTCap).bounding_box(3, "outer");
     
// LCC turn these on?  The bottom is mirrored with the top due to symmetry
// Define the bounding box of the Bottom Toroidal Cap
//     double ksiloB = (*BottomTCap).bounding_box(0, "outer");
//     double etaloB = (*BottomTCap).bounding_box(1, "outer");
//     double ksihiB = (*BottomTCap).bounding_box(2, "outer");
//     double etahiB = (*BottomTCap).bounding_box(3, "outer");
     

// Coil currents at zero flux state:
// These come from ZD "PF_Magnet.doc"
double ccoil1[18] = {1.24,	0.47,	1.85,	7.08,	4.29,	4.05,
					 3.87,	4.20,	4.02,	3.64,	3.02,	2.27,
					 1.35,	0.15,  -1.35,  -3.11,  -4.99,  -6.83};
                     
double ccoil2[18] = {1.03,	1.76,	3.23,	5.85,	3.76,	3.78,
					 3.81,	4.36,	4.37,	4.14,	3.59,	2.86,
					 1.91,	0.57,  -1.17,  -3.32,  -5.74,  -8.15};
                     
double ccoil[18];


// Rate of change of coil current with flux state:                  
double dcdp[18] = {1.40,	0.47,	0.54,	2.44,	0.87,	0.67,
				   0.52,	0.49,	0.31,	0.19,	0.10,	0.02,
				   -0.03,  -0.07,  -0.06,	0.02,	0.14,	0.29};

// Current limiters:
// ... upper limiter of the current
   double qmax; fin_mag.seekg(9*iskip_mag + itxt_mag +1); fin_mag >> qmax;

// ... lower limiter of the current
   double qmin; fin_mag.seekg(10*iskip_mag + itxt_mag +1); fin_mag >> qmin;


// Determine PF coil currents at given q95 and scale to plasma current:
   double qmhd = (*Machine).plasma_edge_safety_factor();
   double aip  = (*Machine).plasma_current();

      for(int i=0; i<18; i++) // 0-17 are all PF coils above midplane
      { if(qmhd >= qmax)
         	{ ccoil[i] = ccoil2[i]*(aip/12.8); }
        if(qmhd <= qmin)
        	{ ccoil[i] = ccoil1[i]*(aip/12.8); }
         if(qmhd > qmin && qmhd < qmax)
         	{ ccoil[i] = ccoil1[i] + ((ccoil2[i]-ccoil1[i])/(qmax-qmin))*(qmax-qmhd);
              ccoil[i] = ccoil[i]*(aip/12.8); } }


// Calculate flux swing required to rampup to Ip
      double r     = (*Machine).plasma_major_radius();
      double asp   = (*Machine).plasma_major_radius() /
                     (*Machine).plasma_minor_radius(); // plasma aspect ratio            
      double ak    = (*Machine).plasma_elongation();   // plasma elongation
      double betap = (*Machine).plasma_poloidal_beta();

      double e   = 1.0/asp;
      double ali = 0.50; // should be in input file...
      double ce  = 0.6;

// Reference ZD PF_Magnet.doc
      double a1 = (1.0 + 1.81*pow(e, 0.5) + 2.05*e)*log(8.0/e)-(2.0 + 9.25*pow(e, 0.5) - 1.21*e);
      double a2 = 0.73*pow(e, 0.5)*(1.0 + 2.0*pow(e, 4.0) - 6.0*pow(e, 5.0) + 3.7*pow(e, 6.0));
      double a3 = 1.0 + 0.98*pow(e, 2.0) + 0.49*pow(e, 4.0) + 1.47*pow(e, 6.0);
      double a4 = 0.25*e*(1.0 + 0.84*e - 1.44*pow(e, 2.0));
      
      double als    = a1*(1.0 - e)/(1.0 - e + a2*ak);
      double am     = pow(1.0 - e, 2.0)/(pow(1.0 - e, 2.0)*a3 + a4*pow(ak, 0.5));
      double alext  = als - (am/4.0)*(log(8.0/(e*pow(ak, 0.5))) + betap + ali/2.0 - 3.0/2.0);
      double alext2 = als;
      double dfcsx  = 1.2566*aip*r*(alext  + (ali/2.0) + ce);
      double dfcsx2 = 1.2566*aip*r*(alext2 + (ali/2.0) + ce);
      

// Calculate forward and back bias coil currents
      double ccoilf[18], ccoilb[18];
      for(int i=0; i<18; i++)
      { ccoilf[i] = ccoil[i]     - (dcdp[i]*dfcsx2/2.0)/25.0;
        ccoilb[i] = ccoil[i]/4.0 + (dcdp[i]*dfcsx2/2.0)/25.0; }


// Determine maximum current in each PF coil
      double ccoilmx[18];
      for(int i=0; i<18; i++)
      { ccoilmx[i] = max( max(ccoilb[i], ccoilf[i]), ccoil[i] ); }


// Calculate coil thickness based on jSC central solenoid (first 4 coils in array)
// outer radius of the central solenoid calculated as inner radius of the
// bucking cylinder minus the gap between central solenoid and bucking cylinder

	 double csbcgap = (*Machine).inboard_generic_build(23); // central solenoid to bucking cyl gap
     double rpfout  = ksiloBC - csbcgap; 

// ... maximum allowed field on PF coils
      double bmaxlimpf; fin_mag.seekg(11*iskip_mag + itxt_mag +1); fin_mag >> bmaxlimpf;
      double ccoilcsmx = max( max(ccoilmx[1], ccoilmx[2]), max(ccoilmx[3], ccoilmx[4]) );
      ccoilcsmx = fabs(ccoilcsmx);
      double bmaxpf = dfcsx/(2.0*pi*pow(rpfout, 2.0));
      
      if(bmaxpf > bmaxlimpf)
      { cout << "bmaxpf = " << bmaxpf << " which is greater than the limit" << bmaxlimpf << "\n";
        (*Machine).get_data_point_filter(3, 3); return; }
      
// maximum allowed current density in superconductor
      double xjsclimpf; 
      if(bmaxpf > 12.0)
      	{ xjsclimpf = -2.1e8*bmaxpf + 3.78e9; } 
      else
      	{ xjsclimpf = -4.0e7*bmaxpf + 1.74e9; }
      double hcs = 2.0*etahiBC; // bucking cylinder height


// Areas: 
   double Asc = ccoilcsmx*1.0e6/xjsclimpf; // total area of the superconductor


// ... maximum stress in protective sheath around strands
   double sigma_sheath; fin_mag.seekg(0*iskip_mag + itxt_mag +1); fin_mag >> sigma_sheath;
   double Asheath = ccoilcsmx*1.0e6*bmaxpf*rpfout/sigma_sheath;  // total area of the sheath


// ... quench protection parameter
   double J2_tau; fin_mag.seekg(12*iskip_mag + itxt_mag +1); fin_mag >> J2_tau;


// ... quench time
   double tau; fin_mag.seekg(13*iskip_mag + itxt_mag +1); fin_mag >> tau;  
   double Acu = ccoilcsmx*1.0e6/pow(2.0*J2_tau/tau,0.5);  // total area of the Cu


// ... volume fraction of He coolant
   double f_He; fin_mag.seekg(3*iskip_mag + itxt_mag +1); fin_mag >> f_He;
   double Ahe = (Asc+Acu)*f_He; // total area of the coolant (He)


// total area of the stainless steel structure
   double Astructure = ccoilcsmx*1.0e6*bmaxpf*rpfout/sigma;  
                                                               

// ... volume fraction of the insulator
   double f_ins; fin_mag.seekg(4*iskip_mag + itxt_mag +1); fin_mag >> f_ins;                                                          
   double Ains = (Asc + Asheath + Acu + Ahe + Astructure)*f_ins/(1.0 - f_ins); 
      
      
// total area of the insulator
   double Atotal = Asc + Acu + Asheath + Astructure + Ahe + Ains;
      
      
// Area fractions
      double fstarea   = (Asheath+Astructure)/Atotal;
      double fcuarea   = Acu/Atotal;
      double fcoolarea = Ahe/Atotal;
      double finsarea  = Ains/Atotal;
      double fscarea   = 1.0 - fstarea - fcuarea - fcoolarea - finsarea;
      double dr        = 8.0*ccoilcsmx*1.0e6/(hcs*fscarea*xjsclimpf);
      double rpfin     = rpfout - dr;
//cout << "rpfout = " << rpfout << "\n"; // ~2.41
//cout << "dr = " << dr << "\n"; // ~0.139
//cout << "ccoilcsmx = " << ccoilcsmx << "\n"; // ~8.01
//cout << "hcs = " << hcs << "\n"; // ~10.04
//cout << "fscarea = " << fscarea << "\n"; // ~0.028
//cout << "xjsclimpf = " << xjsclimpf << "\n"; // ~1.66e9

      	if(rpfin <= 0.0)
      		{ cout << "rpfin in PF coil design is negative, stop\n"; return; }
			// Problem - sometimes this goes negative and fails.  Problem stems from a build
			// issue. Run Matlab loadplot.m to see which component is out of proportion.
			// Can also check Part.cpp ~1010 for list of part IDs.
			// Latest problem was the addition of "string blanketmode" in show()


//  Outer PF coils (coils 5 to 18 in array)
      double dxmx = -1.0e6; // maximum thickness among all the PF coils (initiated)
      double dx[18];
      
      for(int i=0; i<18; i++)
      {  if(i<4){ dx[i] = dr; }
         else
         { dx[i] = pow(fabs(ccoilmx[i]*1.0e6)/(fscarea*xjsclimpf), 0.5);}
         	if(dx[i] > dxmx) dxmx = dx[i]; }


// Define the end-point of the PF coil distribution:
   double x_end_point =  ksihiT;


// Define the PF coil coordinates:    
   double rpf[18], zpf[18];
   string position[18];
   int lttc =  (*TopTCap).contour_size("outer");
         for(int i=0; i<18; i++)
         { if(i<4)
           { rpf[i] = rpfout;
             zpf[i] = etahiBC*i/4.0;
             position[i] = "left"; }
           else
           { double spacing = (ksihiT-ksiloT)/14;

// Locations of the PF coils are determined in the following manner:
// ... all the coils on the central solenoid are on the inner side of the bucking
//     cylinder (i<4);
// ... all the coils on the to top of the toroidal cap are equally spaced. If the
//     slope of the cap surface is going up, they are on the left side of the
//     spacer. If the slope is going down, they are on the right side.
// ... the spacer rpf[i] is placed by dividing the width of the toroidal cap by 14,
//     which is the number of coils on the top. Starting from the 12th coil, the
//     spacer is shifted to the left by dx[i].

           rpf[i] = ksiloT + spacing*(i-3);
           if( i > 9 ) rpf[i] = rpf[i] - dx[i]/2.0;

           double dist = 1.0e6;
           int ipfc;
                for(int j=0; j<lttc; j++)
                {    double xttc = (*TopTCap).contour(j, 0, "outer");
                     double yttc = (*TopTCap).contour(j, 1, "outer");
                           if( dist > fabs(xttc-rpf[i]) &&
                               xttc < (*TopTCap).contour(j+1, 0, "outer") )
                          {    dist = fabs(xttc-rpf[i]);
                               ipfc = j;
                               if ( yttc >= (*TopTCap).contour(j-1, 1, "outer")  )
                               { position[i] = "left"; }
                               else
                               { position[i] = "right"; }
                          } }
               zpf[i] =  (*TopTCap).contour(ipfc, 1, "outer");
           } }


// Generate the contours of the PF coils
for(int k=0; k<18; k++) // loop over top 18 PF coils
{
// ... reference points ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 5;
     double xs[5], ys[5];
     
     
// Central solenoid:
     if( k < 4 )
     { xs[0] = rpf[k] - dx[k];
       ys[0] = zpf[k];

       xs[1] = rpf[k] - dx[k];
       ys[1] = zpf[k] + etahiBC/4;

       xs[2] = rpf[k];
       ys[2] = zpf[k] + etahiBC/4;

       xs[3] = rpf[k];
       ys[3] = zpf[k];

       xs[4] = xs[0]+dx[k]/100.0;
       ys[4] = ys[0]; }
      
      
// PF coil touching the left side of the TF contour  
      if( k >= 4 && position[k] == "left")
    { xs[0] = rpf[k] - dx[k];
      ys[0] = zpf[k];

      xs[1] = rpf[k] - dx[k];
      ys[1] = zpf[k] + dx[k];

      xs[2] = rpf[k];
      ys[2] = zpf[k] + dx[k];

      xs[3] = rpf[k];
      ys[3] = zpf[k];

      xs[4] = xs[0]+dx[k]/100.0;
      ys[4] = ys[0]; }
      
      
// PF coil touching the right side of the TF contour
      if( k >= 4 && position[k] == "right")
    { xs[0] = rpf[k];
      ys[0] = zpf[k];

      xs[1] = rpf[k];
      ys[1] = zpf[k] + dx[k];

      xs[2] = rpf[k] + dx[k];
      ys[2] = zpf[k] + dx[k];

      xs[3] = rpf[k] + dx[k];
      ys[3] = zpf[k];

      xs[4] = xs[0]+dx[k]/100.0;
      ys[4] = ys[0];
      
      if(xs[3] > x_end_point)
      {
//       cout << "PF Coils Too Big \n";
//       (*Machine).get_data_point_filter(2.0); return;
      } }

      for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }     


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }


// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i]     = "line"; }
      connectors_o[nptso] = "none";


// ... Update the coils ...
	double ref_pts_i[100][2], center_i[100][2];
	string connectors_i[100];
	int nptsi = 1;
	int resol_i[100];
	resol_i[0] = 1;

(*Magnet)[k].get_contour_descriptions(nptsi,        nptso,		resol_i,      resol_o,
                                      ref_pts_i,    ref_pts_o,	connectors_i, connectors_o);
double targcurve[10000][2];
int tgl = 0;


// ... Connect the Reference Points into Contours ...
(*Magnet)[k].connect(center_i, center_o, targcurve, tgl);


// ... Establish Contour Depth for volume calculation ...
 double * d_in; double * d_out;
 int lin  = (*Magnet)[k].contour_size("inner");
 int lout = (*Magnet)[k].contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++) d_in[i]   = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Magnet)[k].contour(i, 0, "outer")/16.0;

 (*Magnet)[k].get_depth(d_in, d_out);


// ... Volume and Surface Areas ...
(*Magnet)[k].get_volume();
(*Magnet)[k].get_surface_areas();

delete [] d_in;
delete [] d_out;

} // End of loop over top 18 PF coils


// loop over bottom 18 PF coils
int n = -1;
for(int k=18; k<36; k++) // loop over bottom 18 PF coils
	{ n++; int p = k - 2*n - 1;


// ... ... reference points ... ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 5;
     double xs[5], ys[5];


// Central solenoid:
     if( p < 4 )
    { xs[0] =  rpf[p] - dx[p];
      ys[0] = -zpf[p] - etahiBC/4;

      xs[1] =  rpf[p] - dx[p];
      ys[1] = -zpf[p];

      xs[2] =  rpf[p];
      ys[2] = -zpf[p];

      xs[3] =  rpf[p];
      ys[3] = -zpf[p] - etahiBC/4;

      xs[4] = xs[0]-dx[p]/100.0;
      ys[4] = ys[0]; }


// PF coil touching the left side of the TF contour
      if( p >= 4 && position[p] == "left")
    { xs[0] =  rpf[p] - dx[p];
      ys[0] = -zpf[p] - dx[p];

      xs[1] =  rpf[p] - dx[p];
      ys[1] = -zpf[p];

      xs[2] =  rpf[p];
      ys[2] = -zpf[p];

      xs[3] =  rpf[p];
      ys[3] = -zpf[p] - dx[p];

      xs[4] = xs[0] - dx[p]/100.0;
      ys[4] = ys[0]; }


// PF coil touching the right side of the TF contour
      if( p >= 4 && position[p] == "right")
    { xs[0] =  rpf[p];
      ys[0] = -zpf[p] - dx[p];

      xs[1] =  rpf[p];
      ys[1] = -zpf[p];

      xs[2] =  rpf[p] + dx[p];
      ys[2] = -zpf[p];

      xs[3] =  rpf[p] + dx[p];
      ys[3] = -zpf[p] - dx[p];

      xs[4] = xs[0]-dx[p]/100.0;
      ys[4] = ys[0];
      
      if(xs[3] > x_end_point)
      {
//        cout << "PF Coils Too Big \n";
//       (*Machine).get_data_point_filter(2.0);
//       return;
      } }
      for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }


// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line"; }
      connectors_o[nptso] = "none";


// ... Update the coils ...
	double ref_pts_i[100][2], center_i[100][2];
	string connectors_i[100];
	int nptsi = 1;
	int resol_i[100];
	resol_i[0] = 1;

(*Magnet)[k].get_contour_descriptions(nptsi,        nptso,		resol_i,      resol_o,
                                 	  ref_pts_i,    ref_pts_o,	connectors_i, connectors_o);
double targcurve[10000][2];
int tgl = 0;


// ... Connect the Reference Points into Contours ...
(*Magnet)[k].connect(center_i, center_o, targcurve, tgl);


// ... Establish Contour Depth for volume calculation ...
 double * d_in; double * d_out;
 int lin  = (*Magnet)[k].contour_size("inner");
 int lout = (*Magnet)[k].contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Magnet)[k].contour(i, 0, "outer")/16.0;

 (*Magnet)[k].get_depth(d_in, d_out);


// ... Volume and Surface Areas ...
(*Magnet)[k].get_volume();
(*Magnet)[k].get_surface_areas();

delete [] d_in;
delete [] d_out;

} // end of loop over bottom 18 PF coils

fin_mag.close();

} // eof


// ====================================================================================
// Cryogenic Dome, called "cryodome or cryostat"

void design_cryogenic_dome_shape(DesignPoint * Machine, Part * Cryo, Part * RefBody)
{ 								// RefBody is TFCoil
if( (*Machine).data_point_filter(0) > 0 ) return;
   const double pi = 4.0*atan(1.0);

// Inner Contour of the Cryogenic Dome
   double ksilo = (*RefBody).bounding_box(0, "outer"); // min IB TF coil (x)
   double etalo = (*RefBody).bounding_box(1, "outer"); // bottom TF coil (y)
   double ksihi = (*RefBody).bounding_box(2, "outer"); // max OB TF coil (x)
   double etahi = (*RefBody).bounding_box(3, "outer"); // top of TF coil (y)
   double dwall = (*Machine).vertical_generic_build(21); // cryodome thickness
/* cout << "ksilo = " << ksilo << "\n";
cout << "etalo = " << etalo << "\n";
cout << "ksihi = " << ksihi << "\n";
cout << "etahi = " << etahi << "\n"; */

   int nptsi = 24; // number of points that define contour  
   double xs[30], ys[30];
 
 // Contours reference is Cryogenic_Dome.pdf by Z. Dragojlovic
   xs[0]  = 0.0;
   ys[0]  = 1.5*etahi; // 1.5* OB TF coil dim.
   xs[20] = ksilo + ( ksihi - ksilo )*1.1;
   ys[20] = etahi*1.1;
   
      for(int i=1; i<20; i++)
    { xs[i] = xs[0] + i*xs[20]/20;
      ys[i] = ys[20] + (ys[0] - ys[20])*2.0/pi*atan(2.0*fabs( xs[20] - xs[i])); }
      
      xs[21] = xs[20];
      ys[21] = 1.5*etalo;
      xs[22] = 0.001;
      ys[22] = ys[21];      
      xs[23] = 0.0;
      ys[23] = ys[21];
      
      double ref_pts_i[100][2];     
      for(int i=0; i<nptsi; i++)
     { ref_pts_i[i][0] = xs[i];
       ref_pts_i[i][1] = ys[i]; }
   

// ... ... resolution ... ...
    int resol_i[100];
    for(int i=0; i< nptsi; i++)
    { resol_i[i] = 20; }
      resol_i[nptsi-1] = 0;


// ... ... connectors ... ...
    string connectors_i[100];
    for(int i=0; i<nptsi; i++)
    { connectors_i[i] = "line"; }
      connectors_i[nptsi-1] = "none";
    
    
// Outer Contours of the Cryogenic Dome
   int nptso = 24;
   xs[0] = 0.0;
   ys[0] = 1.5*etahi + dwall;

   xs[20] = ksilo + ( ksihi - ksilo )*1.1 + dwall;
   ys[20] = etahi*1.1;

   for(int i=1; i<20; i++)
    { xs[i] = xs[0] + i*xs[20]/20;
      ys[i] = ys[20] + (ys[0] - ys[20])*2.0/pi*atan(2.0*fabs(xs[20] - xs[i])); }
      
      xs[21] = xs[20];
      ys[21] = 1.5*etalo-dwall;
      xs[22] = 0.001;
      ys[22] = ys[21];
      xs[23] = 0.0;
      ys[23] = ys[21];

      double ref_pts_o[100][2];
      for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }


// ... ... resolution ... ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }
    resol_o[nptso-1] = 0;


// ... ... connectors ... ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line"; }
    connectors_o[nptso-1] = "none";


// ... Update the Cryogenic Dome ...
(*Cryo).get_contour_descriptions(nptsi, 		nptso,		  resol_i, 
								 resol_o,   	ref_pts_i, 	  ref_pts_o,
                                 connectors_i,  connectors_o);


// Define the Inner Target Curve For Blanket Shape
double targcurve[10000][2];
int tgl = 0;
    

// ... Connect the Reference Points into Contours ...
double center_i[100][2], center_o[100][2];
(*Cryo).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 double * d_in; double * d_out;
 int lin =  (*Cryo).contour_size("inner");;
 int lout = (*Cryo).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];


// Main body of the cryodome
 for(int i=0; i<lin; i++)   d_in[i] = 2.0*pi*(*Cryo).contour(i, 0, "inner")/16.0;
 for(int i=0; i<lout; i++) d_out[i] = 2.0*pi*(*Cryo).contour(i, 0, "outer")/16.0;
	(*Cryo).get_depth(d_in, d_out);

// Volume and Surface Areas
	(*Cryo).get_volume();
	(*Cryo).get_surface_areas();	
//cout << "Cryostat volume (G.cpp) = " << (*Cryo).volume() << "\n"; // unit vol ~65 m3

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Maintenance Ports

void design_maintenance_port(DesignPoint *Machine, Part *MPort, 
							 Part (*Vessel)[4], Part *Cryo)
{ 
if( (*Machine).data_point_filter(0) > 0 ) return;

// ... Define the Inner Contour ...
     double ksihi = (*Cryo).bounding_box(2, "outer"); // OB cryodome (x)
     double OBvvx = (*Vessel)[0].bounding_box(2, "outer"); // max OB TF coil (x)

// ... ... reference points ... ...
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];


// ... Define the Outer Contour ...
// ... ... reference points ... ...
     double ref_pts_o[100][2], center_o[100][2];
     int nptso = 4;
     double xs[10], ys[10];

// Reference points:
// first: ref start or end; 0=start ref pt, then 0(+1)=end
// second: x=[0], y=[1]
// third: choose inner or outer side of component
if (blanketmode == "dcll") {
	 OBvvx = (*Vessel)[1].bounding_box(2, "outer"); // max OB TF coil (x)
	 
     xs[0] = (*Vessel)[1].reference_point(3, 0, "outer");
     ys[0] = (*Vessel)[1].reference_point(3, 1, "outer");

     xs[1] = (*Vessel)[1].reference_point(2, 0, "outer");
     ys[1] = (*Vessel)[1].reference_point(2, 1, "outer");

     xs[2] = ksihi;
     ys[2] = (*Vessel)[1].reference_point(2, 1, "outer");

     xs[3] = ksihi;
     ys[3] = (*Vessel)[1].reference_point(3, 1, "outer"); }

else { 
     xs[0] = (*Vessel)[0].reference_point(3, 0, "outer");
     ys[0] = (*Vessel)[0].reference_point(3, 1, "outer");

     xs[1] = (*Vessel)[0].reference_point(2, 0, "outer");
     ys[1] = (*Vessel)[0].reference_point(2, 1, "outer");

     xs[2] = ksihi;
     ys[2] = (*Vessel)[0].reference_point(2, 1, "outer");

     xs[3] = ksihi;
     ys[3] = (*Vessel)[0].reference_point(3, 1, "outer");     
     }
     
/*cout << "xs[0] = " << xs[0] << "\n";
cout << "ys[0] = " << ys[0] << "\n";
cout << "xs[1] = " << xs[1] << "\n";
cout << "ys[1] = " << ys[1] << "\n";
cout << "xs[2] = " << xs[2] << "\n";
cout << "ys[2] = " << ys[2] << "\n";
cout << "xs[3] = " << xs[3] << "\n";
cout << "ys[3] = " << ys[3] << "\n";*/

   
     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

/* ... Close the loop by adding a symmetric profile in the clockwise direction: ... */
      for(int i=0;  i<nptso; i++)
      { ref_pts_o[nptso+i][0] =  ref_pts_o[nptso-i-2][0];
        ref_pts_o[nptso+i][1] = -ref_pts_o[nptso-i-2][1]; }
        nptso = 2*nptso-2;

// ... resolution ...
    int resol_o[100];
    for(int i=0; i< nptso; i++)
    { resol_o[i] = 20; }

// ... connectors ...
    string connectors_o[100];
    for(int i=0; i<nptso; i++)
    { connectors_o[i] = "line";
      if(i==0 || i==5) connectors_o[i] = "target"; }


// Update the Maintenance Port
(*MPort).get_contour_descriptions(nptsi, nptso, resol_i, resol_o,
                                  ref_pts_i, ref_pts_o, connectors_i, connectors_o);

// Define the Target Curve For the Maintenance Port along OB VV
double targcurve[10000][2];  int tgl;
if (blanketmode == "dcll") { 
	   tgl = (*Vessel)[1].contour_size("outer");
	for(int i=0; i<tgl; i++)
  		{ targcurve[i][0] = (*Vessel)[1].contour(i, 0, "outer");
    	  targcurve[i][1] = (*Vessel)[1].contour(i, 1, "outer"); } }
else { tgl = (*Vessel)[0].contour_size("outer"); // =450, # pts along OB VV
	for(int i=0; i<tgl; i++)
		{ targcurve[i][0] = (*Vessel)[0].contour(i, 0, "outer");
    	  targcurve[i][1] = (*Vessel)[0].contour(i, 1, "outer"); } }


// Connect the Reference Points into Contours
(*MPort).connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 const double pi = 4.0*atan(1.0);
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*MPort).contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];
 
// Duct wall thickness seems to be VV thickness/4, reference?
 	double ductwall = (*Machine).vertical_generic_build(11)/4.0; // =0.1, 1/4th VV thickness
 	for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 	for(int i=0; i<lout; i++) d_out[i] = ductwall;
    (*MPort).get_depth(d_in, d_out);


// Volume and Surface Areas
// ... volume of one side wall of the maintenance port
     (*MPort).get_volume();
     double vSideWall = (*MPort).volume();

// ... volume of the top and bottom walls of the maintenance port
     double rTop     = (*MPort).bounding_box(0, "outer"); // top-right (x) of MP+VV intersection
     double lTop     = ksihi - rTop; // length of MPort (top)
     double vTopWall = 2.0*rTop*sin(pi/16.0)*lTop*ductwall; // vol top wall
	 // straight sides, not a wedge

// ... volume of the outer door of the maintenance port
     double etaMPhi    = (*MPort).bounding_box(3, "outer"); // top-right (y) of MP+VV intersection
     double etaMPlo    = (*MPort).bounding_box(1, "outer"); // bottom-right (y) of MP+VV intersection
     double vOuterDoor = 2.0*rTop*sin(pi/16.0)*(etaMPhi-etaMPlo)*ductwall; 
	 // 2.0*sin(pi/16.0)*R is arclength at R 
	 
// ... amount of volume to add to already included side wall
     double vPlus = 2.0*vTopWall + vSideWall + vOuterDoor;
     (*MPort).add_to_volume(vPlus); // adds vPlus to *MPort part

// ... Volume enclosed w/in MPort
	 double MPortArea = 2.0*( 0.5*(OBvvx-rTop) + (ksihi-OBvvx)*etaMPhi ); // approx area of MPort
	 double Rrevolve  = (OBvvx-rTop)*0.707+rTop;
	 double MPortVol  = MPortArea*2.0*sin(pi/16.0)*Rrevolve; // 1/16th; 0.707 ensures approx is ok for triangle area
     (*MPort).get_enclosed_volume(MPortVol); // adds MPortVol to function to be called later

//cout << "MPortVol = " << 16.0*(*MPort).enclosed_volume() << "\n";
//cout << "MPortVol = " << 16.0*MPortVol << "\n"; // ~2320 m3
//cout << "(*MPort).volume() = " << (*MPort).volume() << "\n"; // total vol 1/16th (not enclosed vol)
//cout << "ductwall = " << ductwall << "\n";
//cout << "etaMPlo = " << etaMPlo << "\n";
//cout << "vOuterDoor = " << vOuterDoor << "\n";

delete [] d_in;
delete [] d_out;


} // eof


// ====================================================================================
// Upper Vacuum Duct, Bottom Wall
// Needs modification - does not scale to machine

void design_upper_vv_duct_bottom_wall(DesignPoint * Machine, 
	Part (*Duct)[3], Part (*Vessel)[4], Part * Magnet, Part * Cryo)	
{ 
if( (*Machine).data_point_filter(0) > 0 ) return; 
   const double pi = 4.0*atan(1.0);

// Relevant bounding box coordinates of the vacuum vessel
	double etahiVV, ksihiVV, ksiloVV;
if (blanketmode == "dcll") { 
	   etahiVV = (*Vessel)[2].bounding_box(3, "outer");
	   ksihiVV = (*Vessel)[2].bounding_box(2, "outer");
	   ksiloVV = (*Vessel)[2].bounding_box(0, "outer"); }
else { etahiVV = (*Vessel)[0].bounding_box(3, "outer");   // top of vert VV (y)
	   ksihiVV = (*Vessel)[0].bounding_box(2, "outer");   // OB VV (x)
	   ksiloVV = (*Vessel)[0].bounding_box(0, "outer"); } // IB VV (x)
double ksimedVV = 0.4*(ksiloVV + ksihiVV); // middle of VV
/*cout << "etahiVV = " << etahiVV << "\n";
cout << "ksihiVV = " << ksihiVV << "\n";
cout << "ksiloVV = " << ksiloVV << "\n"; */

// Relevant bounding box coordinates of the outer cylindrical surface of the cryostat
   double ksihiCryo = (*Cryo).bounding_box(2, "outer");

// Ordinate of the highest point of the inner surface of the TF magnet
   double etahiMagnet = (*Magnet).bounding_box(3, "inner");

// Largest vertical gap between the vacuum vessel and the TF coil
   double hVVTF = etahiMagnet - etahiVV;

// Width of the duct in the vertical plane
   double ductwidth = 0.15*(ksihiVV - ksiloVV); // 0.25 reference?


// Inner contour descriptors are declared but not filled with any data since
// the inner contour does not exist in this case
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];

// Define the top surface of the vacuum duct by using the procedure for the 
// outer contour of the generic "Part" entity.
 double ref_pts_o[100][2], center_o[100][2];
 int nptso = 5; // number of reference points
 double xs[5], ys[5];
 
      xs[0] = ksimedVV + ductwidth;
      ys[0] = etahiVV;

      xs[1] = ksimedVV + ductwidth;
      ys[1] = etahiVV + 0.25*hVVTF;

      ys[2] = etahiVV + 2.0*hVVTF - ductwidth;
      xs[2] = xs[1] + (ys[2]-ys[1])/tan(pi/4.0);

      xs[3] = (2.0/3.0)*ksihiCryo;
      ys[3] = ys[2];

      xs[4] = 1.1*ksihiCryo; // horiz distance beyond cryostat, was 2.0*
      ys[4] = ys[2];

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i]; }

// ... ... resolution ... ...
     int resol_o[100];
     resol_o[0] = 10;
     resol_o[1] = 20;
     resol_o[2] = 50;
     resol_o[3] = 150;
     resol_o[4] = 0;

// ... ... connectors ... ...
    string connectors_o[100];
    connectors_o[0] = "line";
    connectors_o[1] = "line";
    connectors_o[2] = "line";
    connectors_o[3] = "line";
    connectors_o[4] = "none";


// Update the Upper Duct Wall Contour
(*Duct)[0].get_contour_descriptions(nptsi, nptso, resol_i, resol_o,
                                    ref_pts_i, ref_pts_o, connectors_i, connectors_o);

// Define the Inner Target Curve For the Upper Duct Wall
double targcurve[10000][2]; int tgl=0;

// Connect the Reference Points into Contours
(*Duct)[0].connect(center_i, center_o, targcurve, tgl);

// Establish Contour Depth for volume calculation ...
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*Duct)[0].contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = ductwidth/2.0;

 (*Duct)[0].get_depth(d_in, d_out);


// Volume and Surface Areas
// Duct wall thickness seems to be VV thickness/4.0, reference?
	double ductwall = (*Machine).vertical_generic_build(11)/4.0; // 1/4th VV thickness
	(*Duct)[0].get_surface_areas();
	(*Duct)[0].get_shell_volume(ductwall);
	(*Duct)[0].get_cross_section(0, ductwall, ductwidth/2.0);
//cout << "duct wall thickness = " << ductwall << "\n";

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Upper Vacuum Duct, Top Wall

void design_upper_vv_duct_top_wall(DesignPoint * Machine, 
	Part (*Duct)[3], Part (*Vessel)[4], Part * Magnet, Part * Cryo)	
{ 
if( (*Machine).data_point_filter(0) > 0 ) return;
	const double pi = 4.0*atan(1.0);

// Relevant bounding box coordinates of the vacuum vessel
	double etahiVV, ksihiVV, ksiloVV;
if (blanketmode == "dcll") { 
	   etahiVV = (*Vessel)[2].bounding_box(3, "outer");
	   ksihiVV = (*Vessel)[2].bounding_box(2, "outer");
	   ksiloVV = (*Vessel)[2].bounding_box(0, "outer"); }
else { etahiVV = (*Vessel)[0].bounding_box(3, "outer");
	   ksihiVV = (*Vessel)[0].bounding_box(2, "outer");
	   ksiloVV = (*Vessel)[0].bounding_box(0, "outer"); }
double ksimedVV = 0.4*(ksiloVV + ksihiVV); // middle of VV


// Relevant bounding box coordinates of the outer cylindrical surface of the cryostat
 double ksihiCryo = (*Cryo).bounding_box(2, "outer");

// Ordinate of the highest point of the inner surface of the TF magnet
 double etahiMagnet = (*Magnet).bounding_box(3, "inner");

// Largest vertical gap between the vacuum vessel and the TF coil
 double hVVTF = etahiMagnet - etahiVV;

// Duct total width in the vertical plane
 double ductwidth = 0.15*(ksihiVV - ksiloVV); // 0.25 ref?

// Inner contour descriptors are declared but not filled with any data since
// the inner contour does not exist in this case
     double ref_pts_i[100][2], center_i[100][2];
     int nptsi=1;
     int resol_i[100];
     string connectors_i[100];


// Define the top surface of the vacuum duct by using the procedure for the 
// outer contour of the generic "Part" entity.
 double ref_pts_o[100][2], center_o[100][2];
 int nptso = 3; // number of reference points
 double xs[3], ys[3];
 
      xs[0] = ksimedVV;
      ys[0] = etahiVV;

      xs[1] = (2.0/3.0)*ksihiCryo; // 2/3?
      ys[1] = etahiVV + 2.0*hVVTF;

      xs[2] = 1.1*ksihiCryo; // horiz distance beyond cryostat, was 2.0*
      ys[2] = etahiVV + 2.0*hVVTF;

     for(int i=0; i<nptso; i++)
     { ref_pts_o[i][0] = xs[i];
       ref_pts_o[i][1] = ys[i];}

// ... ... resolution ... ...
     int resol_o[100];
     resol_o[0] = 200;
     resol_o[1] = 200;
     resol_o[2] = 0;
    
// ... ... connectors ... ...
    string connectors_o[100];
    connectors_o[0] = "target";
    connectors_o[1] = "line";
    connectors_o[2] = "none";


// Update the Upper Duct Wall Contour
(*Duct)[1].get_contour_descriptions(nptsi, nptso, resol_i, resol_o,
                                    ref_pts_i, ref_pts_o, connectors_i, connectors_o);

// Define the Inner Target Curve For the Upper Duct Wall
double targcurve[10000][2];
int tgl = (*Duct)[0].contour_size("outer");
  for(int i=0; i<tgl; i++)
  { targcurve[i][0] = (*Duct)[0].contour(i, 0, "outer");
    targcurve[i][1] = (*Duct)[0].contour(i, 1, "outer"); }

// Connect the Reference Points into Contours
(*Duct)[1].connect(center_i, center_o, targcurve, tgl);


// Establish Contour Depth for volume calculation ...
 double * d_in; double * d_out;
 int lin = 1;
 int lout = (*Duct)[1].contour_size("outer");

 d_in  = new double [lin];
 d_out = new double [lout];

 for(int i=0; i<lin; i++)   d_in[i] = 0.0;
 for(int i=0; i<lout; i++) d_out[i] = ductwidth/4.0;

 (*Duct)[1].get_depth(d_in, d_out);


// Volume and Surface Areas
	double ductwall = (*Machine).vertical_generic_build(11)/4.0; // 1/4th VV thickness
	(*Duct)[1].get_surface_areas();
	(*Duct)[1].get_shell_volume(ductwall);
	(*Duct)[1].get_cross_section(0, ductwall, ductwidth/4.0);

delete [] d_in;
delete [] d_out;

} // eof


// ====================================================================================
// Upper Vacuum Duct, Side Wall
// LCC - I don't think this function creates any contours for the side wall...
// Plus, this is giving problems in DP.cpp

void design_upper_vv_duct_side_wall(DesignPoint * Machine, Part (*Duct)[3])
{
// Vacuum Duct has a trapezoidal cross section.
// Top side length is denoted as "a", bottom side as "b", lateral sides as "c"
// and height as "h". The thickness of the vacuum duct wall is denoted as "d".

// Sides of the Trapezoid
double a = (*Duct)[1].cross_section_depth(0); // top side of the trapezoidal cross section
double b = (*Duct)[0].cross_section_depth(0); // bottom side of the trapezoidal cross section

double etahiTopWall    = (*Duct)[1].bounding_box(3, "outer");
double etahiBottomWall = (*Duct)[0].bounding_box(3, "outer");

double h = (etahiTopWall - etahiBottomWall); 			  // height of the trapezoid
double c = pow( pow(h, 2.0) + pow( (b-a)/2.0, 2.0), 0.5); // lateral side of the trapezoidal cross section

// Wall Thickness
double dwall = (*Duct)[0].cross_section_width(0);

// Total Lengths of the Top and Bottom Wall of the Vacuum Duct
double lTW = (*Duct)[1].volume()/(dwall*a);
double lBW = (*Duct)[0].volume()/(dwall*b);

// Volume of the Side Wall of the Vacuum Duct
double vSW = 0.5*(lTW + lBW)*c*dwall;
(*Duct)[2].add_to_volume(vSW); // only place Duct[2] is used here

} // eof