// Member Functions of the Class "Part"
// Blanket independent
// LCC cleaned up formatting 9/15/10, could condense more...

#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <fstream>
#include <cmath>
#include <string>

using std::string;
using namespace std;

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


// ====================================================================================
// Contour Descriptions by Reference Points and Connectors
void Part::get_contour_descriptions(int ptsi,            int ptso,
                                    int resi[100],       int reso[100],
                                    double refi[100][2], double refo[100][2],
                                    string coni[100],    string cono[100])
{  num_ref_pts_inner = ptsi;
   num_ref_pts_outer = ptso;
   
    for(int i=0; i<ptsi; i++)
    { connectors_inner[i] = coni[i];
      resolution_inner[i] = resi[i];
      for(int j=0; j<2; j++)
      { ref_pts_inner[i][j] = refi[i][j]; } }
    
    for(int i=0; i<ptso; i++)
    { connectors_outer[i] = cono[i];
      resolution_outer[i] = reso[i];
      for(int j=0; j<2; j++)
      { ref_pts_outer[i][j] = refo[i][j]; } }
} // eof


// ====================================================================================
// Return Bounding Box of a Contour
// ... returns xlo, ylo, xhi, yhi for bx = 0,1,2,3
double Part::bounding_box(int bx, string flag)
{
      double xlo =  1.0e28;
      double xhi = -1.0e28;
      double ylo =  1.0e28;
      double yhi = -1.0e28;

      if(flag=="inner")
      { for(int i=0; i<len_ctr_inner; i++)
       { xlo = (contour_inner[i][0] < xlo) ? contour_inner[i][0] : xlo;
		 xhi = (contour_inner[i][0] > xhi) ? contour_inner[i][0] : xhi;
		 ylo = (contour_inner[i][1] < ylo) ? contour_inner[i][1] : ylo;
		 yhi = (contour_inner[i][1] > yhi) ? contour_inner[i][1] : yhi; } }
	   else
      { for(int i=0; i<len_ctr_outer; i++)
       { xlo = (contour_outer[i][0] < xlo) ? contour_outer[i][0] : xlo;
	     xhi = (contour_outer[i][0] > xhi) ? contour_outer[i][0] : xhi;
		 ylo = (contour_outer[i][1] < ylo) ? contour_outer[i][1] : ylo;
		 yhi = (contour_outer[i][1] > yhi) ? contour_outer[i][1] : yhi; } }

	if(bx==0) return xlo;
	if(bx==1) return ylo;
	if(bx==2) return xhi;
	if(bx==3) return yhi;

    return 0.0;
} // eof


// ====================================================================================
// Return Bounding Box of a Contour at Equatorial Plane (for TETRA format radial build)
// ... returns xlo for bx = 0 and xhi for bx= 1
double Part::bounding_box_equatorial(int bx, string flag)
{
      double xlo =  1.0e28;
      double xhi = -1.0e28;

      if(flag=="inner")
      { for(int i=0; i<len_ctr_inner; i++)
       { if(contour_inner[i][1] >= -1.e-6 & contour_inner[i][1] <= 1.e-6)
        { xlo = (contour_inner[i][0] < xlo) ? contour_inner[i][0] : xlo;
		  xhi = (contour_inner[i][0] > xhi) ? contour_inner[i][0] : xhi; } } }
	  else 
	  { for(int i=0; i<len_ctr_outer; i++)
   	   { if(contour_outer[i][1] >= -1.e-6 & contour_outer[i][1] <= 1.e-6)
        { xlo = (contour_outer[i][0] < xlo) ? contour_outer[i][0] : xlo;
		  xhi = (contour_outer[i][0] > xhi) ? contour_outer[i][0] : xhi; } } }

	if(bx==0) return xlo; // IB?
	if(bx==1) return xhi; // OB?
	
	return 0.0; } // eof


// ====================================================================================
// Return A Contour
double Part::contour( int num, int coord, string flag)
	{ if(flag == "inner") return contour_inner[num][coord];
      if(flag == "outer") return contour_outer[num][coord];
      return 0.0; } // eof


// ====================================================================================
// Return A Reference Point
double Part::reference_point( int num, int coord, string flag)
	{ if(flag == "inner") return ref_pts_inner[num][coord];
   	  if(flag == "outer") return ref_pts_outer[num][coord];
      return 0.0; } // eof


// ====================================================================================
// Return A Contour Size
int Part::contour_size(string flag)
	{ if(flag == "inner") return len_ctr_inner;
      if(flag == "outer") return len_ctr_outer;
      return 0; } // eof


// ====================================================================================
// Connector lines
void Part::connector_line(int ref_start, int con_start, int con_end, string con_id)
{
// ... coordinates of the starting and ending point ...
// 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
        double xstart, xend, ystart, yend;
        if( con_id == "inner")
        {  xstart =  ref_pts_inner[ref_start][0];
           ystart =  ref_pts_inner[ref_start][1];
           xend =    ref_pts_inner[ref_start+1][0];
           yend =    ref_pts_inner[ref_start+1][1]; } 
        else
        {  xstart =  ref_pts_outer[ref_start][0];
           ystart =  ref_pts_outer[ref_start][1];
           xend =    ref_pts_outer[ref_start+1][0];
           yend =    ref_pts_outer[ref_start+1][1]; }
        
// ... line segment between the points ...
      double length = sqrt( pow( (xend-xstart),2) + pow( (yend-ystart),2) );
      double delta = length/(float(con_end) - float(con_start));
      
      if(con_id == "inner")
      { contour_inner[con_start][0] = xstart;
		contour_inner[con_start][1] = ystart;
         for(int i=0; i<con_end-con_start; i++)
         { contour_inner[con_start+i][0] = xstart + (xend-xstart)/length*delta*i;
           contour_inner[con_start+i][1] = ystart + (yend-ystart)/length*delta*i; } }
      else 
      {	contour_outer[con_start][0] = xstart;
		contour_outer[con_start][1] = ystart;
         for(int i=0; i<con_end-con_start; i++)
         { contour_outer[con_start+i][0] = xstart + (xend-xstart)/length*delta*i;
           contour_outer[con_start+i][1] = ystart + (yend-ystart)/length*delta*i; } }
} // eof


// ====================================================================================
// Connector circles
void Part::connector_circle(double xcen,   double ycen,
                            int ref_start, int con_start,
                            int con_end,   string con_id)
{
// ... coordinates of the starting and ending point ...
        double xstart, xend, ystart, yend;
        if( con_id == "inner")
        {  xstart =  ref_pts_inner[ref_start][0];
           ystart =  ref_pts_inner[ref_start][1];
           xend =    ref_pts_inner[ref_start+1][0];
           yend =    ref_pts_inner[ref_start+1][1]; } 
        else
        {  xstart =  ref_pts_outer[ref_start][0];
           ystart =  ref_pts_outer[ref_start][1];
           xend =    ref_pts_outer[ref_start+1][0];
           yend =    ref_pts_outer[ref_start+1][1]; }
       double radius = sqrt( pow( (xcen-xstart),2) + pow( (ycen-ystart),2) );
       double dx = (xend-xstart)/(float(con_end - con_start));
       double d_plus = ycen - ystart + sqrt( pow(radius,2) - pow( (xstart - xcen),2) );
      

 /* ... decide whether to connect upper or lower part of the circle ... */
      int sgn;     
      if(d_plus == 0.0)
      	{ sgn = 1; }
      else
      	{ sgn = -1; }

     if(con_id == "inner")
     	{ for(int i=0; i<=con_end-con_start; i++)
      	{ contour_inner[con_start+i][0] = xstart + dx*float(i);
          contour_inner[con_start+i][1] = ycen +
      	  float(sgn)*sqrt( pow(radius,2) - pow( (xstart + dx*float(i) - xcen),2)); } }

     if(con_id == "outer")
     	{ for(int i=0; i<=con_end-con_start; i++)
        { contour_outer[con_start+i][0] = xstart + dx*float(i);
          contour_outer[con_start+i][1] =  ycen +
          float(sgn)*sqrt( pow(radius,2) - pow( (xstart + dx*float(i) - xcen),2)); } }
} // eof


// ====================================================================================
// Connectors
void Part::connector_target_curve(double target[10000][2], int targlen,
                               	  int ref_start,           int con_start,
                                  int *con_end,            string con_id)
{
// ... coordinates of the starting and ending point ...
        double xstart, xend, ystart, yend;
        if( con_id == "inner")
        {  xstart =  ref_pts_inner[ref_start][0];
           ystart =  ref_pts_inner[ref_start][1];
           xend =    ref_pts_inner[ref_start+1][0];
           yend =    ref_pts_inner[ref_start+1][1];
        } else
        {  xstart =  ref_pts_outer[ref_start][0];
           ystart =  ref_pts_outer[ref_start][1];
           xend =    ref_pts_outer[ref_start+1][0];
           yend =    ref_pts_outer[ref_start+1][1]; }

        
/* ... get the starting and ending distances from the target curve ... */
    double dstart = 1.0e38;
	double dend   = 1.0e38;
	double distref;
	int istart, iend;

	for(int i=0; i<targlen; i++)
		{ distref = sqrt(pow((xstart-target[i][0]),2) +
                         pow((ystart-target[i][1]),2));
           if(dstart >= distref )
           { dstart = distref;
	   		 istart = i; }
           distref =  sqrt(pow((xend-target[i][0]),2) +
                           pow((yend-target[i][1]),2));            
          if(dend >= distref)
          { dend = distref;
            iend = i; } }
    int imin = (istart < iend) ? istart : iend;
	int imax = (istart > iend) ? istart : iend;
	

/* ... calculate the total length of the segment ... */
    double length = 0.0;
	for(int i = imin; i<=imax; i++)
        { length = length + sqrt(pow((target[i+1][0]-target[i][0]),2) +
                                 pow((target[i+1][1]-target[i][1]),2)); }

/* ... get the designed curve ... */
    double t = 0.0;
	int N = con_start;
       for(int i=imin; i<=imax; i++)
       	{ double tagm;  
		if(i>0)
        	{ tagm = sqrt(pow((target[i+1][0]-target[i-1][0]),2) +
                          pow((target[i+1][1]-target[i-1][1]),2)); }
        else
        { tagm = 2.0*sqrt(pow((target[i+1][0]-target[i][0]),2) +
                          pow((target[i+1][1]-target[i][1]),2)); }

	double nx = - (target[i+1][1]-target[i-1][1])/tagm;
	double ny =   (target[i+1][0]-target[i-1][0])/tagm;
    double dvar;
        
        if(imin == istart)
        	{ dvar = dstart*(length-t)/length + dend*t/length; }
        else
        	{ dvar = dend*(length-t)/length + dstart*t/length; }

        if(con_id=="inner")
        	{ contour_inner[N][0] = target[i][0] + nx*dvar;
		      contour_inner[N][1] = target[i][1] + ny*dvar;	}
		else
        	{ contour_outer[N][0] = target[i][0] + nx*dvar;
		      contour_outer[N][1] = target[i][1] + ny*dvar; }
         N++;
         t = t + tagm/2.0; }
         *con_end = N-1;

      if(imin==iend)
      	{ double swap;
      for(int i=con_start; i<= con_start + (*con_end - con_start)/2; i++)
      	{ if(con_id=="inner")
        	{ swap = contour_inner[i][0];
			  contour_inner[i][0] = contour_inner[*con_end + con_start - i][0];
			  contour_inner[*con_end + con_start - i][0] = swap;
			  
			  swap = contour_inner[i][1];
			  contour_inner[i][1] = contour_inner[*con_end + con_start - i][1];
			  contour_inner[*con_end + con_start - i][1] = swap; }
		  else
        	{ swap = contour_outer[i][0];
			  contour_outer[i][0] = contour_outer[*con_end + con_start -i][0];
			  contour_outer[*con_end + con_start - i][0] = swap;

			  swap = contour_outer[i][1];
			  contour_outer[i][1] = contour_outer[*con_end + con_start - i][1];
			  contour_outer[*con_end + con_start - i][1] = swap; }
       } // end of for loop (i)
       } // end of imin==iend condition


// ... fix the end-points ...
	if(con_id=="inner")
    	{ contour_inner[con_start][0] = xstart;
		  contour_inner[con_start][1] = ystart;
		  contour_inner[*con_end][0] = xend;
		  contour_inner[*con_end][1] = yend; }
	else
        { contour_outer[con_start][0] = xstart;
		  contour_outer[con_start][1] = ystart;
		  contour_outer[*con_end][0] = xend;
		  contour_outer[*con_end][1] = yend; }

/* ... iron out the contour by removing any internal loops ... */
        disentangle(con_end, con_id);
} // eof


// ====================================================================================
// Disentangle the Loops in Contour Caused by Sharp Corners in Target Curve

// translate sharp corners into sharp corners instead of sharp corners into loops
void Part::disentangle(int * n, string cflag)
{
/* ... tag the points that lie on the enclosed loops ... */
	string tag[10000];
    for(int i=0; i<*n-1; i++) tag[i] = "off";
    for(int i=0; i<*n-1; i++)
    	{ double xstart, ystart, xend, yend;
    if(cflag == "inner")
    	{ xstart = contour_inner[i][0];
	  	  ystart = contour_inner[i][1];
	  	  xend = contour_inner[i+1][0];
	  	  yend = contour_inner[i+1][1]; }
	else
		{ xstart = contour_outer[i][0];
      	  ystart = contour_outer[i][1];
	  	  xend = contour_outer[i+1][0];
	  	  yend = contour_outer[i+1][1]; }

	double ksweep =i;
	for(int j=i; j<*n-1; j++)
      { double xs, ys, xe, ye;
      if(cflag == "inner")
      	{ xs = contour_inner[j][0];
	      ys = contour_inner[j][1];
	      xe = contour_inner[j+1][0];
          ye = contour_inner[j+1][1]; }
      else
      	{ xs = contour_outer[j][0];
	      ys = contour_outer[j][1];
	      xe = contour_outer[j+1][0];
          ye = contour_outer[j+1][1]; }
       
	double bigreal = 1.e38;        
    double xint = bigreal;
    double yint = bigreal;
    intersect(xstart, ystart, xend, yend, xs, ys, xe, ye, &xint, &yint, bigreal);

   if(xint != bigreal && yint != bigreal && j > i+1 &&
      xint >= ((xs < xe) 		? xs : xe) 			&& xint <= ((xs > xe) 		? xs : xe) &&
      yint >= ((ys < ye) 		? ys : ye) 			&& yint <= ((ys > ye) 		? ys : ye) &&
      xint >= ((xstart < xend) 	? xstart : xend) 	&& xint <= ((xstart > xend) ? xstart : xend) &&
      yint >= ((ystart < yend) 	? ystart : yend) 	&& yint <= ((ystart > yend) ? ystart : yend))
    	{ ksweep = j+1; break; }
        } // end for loop (j)
        
    for(int j=i+1; j<ksweep; j++) tag[j] = "on";        
        }  // end for loop (i)


/* ... eliminate the points that are tagged ...*/
      int j = -1;
      for(int i=0; i<*n; i++)
      { if(tag[i] == "off" || i == 0 || i == *n-1)
      	{ j++;
        	if(cflag == "inner")
             { contour_inner[j][0] = contour_inner[i][0];
               contour_inner[j][1] = contour_inner[i][1];} 
            else
             { contour_outer[j][0] = contour_outer[i][0];
               contour_outer[j][1] = contour_outer[i][1]; } } }
      *n = j; }  // eof


// ====================================================================================
// Intersection Between Two Vectors
int Part::intersect( double x1s,    double y1s,    double x1e, double y1e,
                     double x2s,    double y2s,    double x2e, double y2e,
                     double * xint, double * yint, double bigreal)
{
      string flag[3];
      double a1, a2, b1, b2;
      
      if(x1e != x1s)
      	{ flag[0] = "on";
          a1 = (y1e-y1s)/(x1e-x1s);
          b1 = y1s - x1s*a1; } 
      else
	  	{ flag[0] = "off";
	  	  a1 = bigreal;
	      b1 = bigreal; }

      if(x2e != x2s)
      	{ flag[1] = "on";
          a2 = (y2e-y2s)/(x2e-x2s);
          b2 = y2s - x2s*a2; }
       else
      	{ flag[1] = "off";
		  a2 = bigreal;
		  b2 = bigreal; }

	  if(a1 != a2)
		{ flag[2] = "on"; }
		else
		{ flag[2] = "off"; }

	  if(flag[0] == "on" && flag[1] == "on" && flag[2] == "on")
		{ *xint = -(b1-b2)/(a1-a2);
	      *yint = a1*(*xint) + b1;
	      return 0; }

	  if(flag[0] == "off" && flag[1] == "on" && flag[2] == "on")
		{ *xint = x1s;
	      *yint = a2*(*xint) + b2;
	      return 0; }

	  if(flag[0] == "on" && flag[1] == "off" && flag[2] == "on")
		{ *xint = x2s;
	      *yint = a1*(*xint) + b1;
	      return 0; }

	  if( (flag[0] == "off" && flag[1] == "off") || flag[2] == "off" )
		{ *xint = bigreal;
          *yint = bigreal;
          return 0; }
          return 0; } // eof


// ====================================================================================
// Connect the Contours
// connects the reference points into a single contour by using the connector functions

void Part::connect(double center_inner[100][2], double center_outer[100][2],
                   double target[10000][2],     int targlen)
{
    int con_start = 0;
    int con_end   = 0;
    len_ctr_inner = 0;
    if(num_ref_pts_inner > 1)
    	{ for(int i=0; i<num_ref_pts_inner; i++)
     	{ if(connectors_inner[i] == "target")
       	{ len_ctr_inner = len_ctr_inner + targlen; } 
    else
        { len_ctr_inner = len_ctr_inner + resolution_inner[i]; } } } 
    else
     	{ len_ctr_inner = 1;}
    
    len_ctr_outer = 0;
    for(int i=0; i<num_ref_pts_outer; i++)
    	{ if(connectors_outer[i] == "target")
        { len_ctr_outer = len_ctr_outer + targlen; } 
    else
        { len_ctr_outer = len_ctr_outer + resolution_outer[i]; } }
        
    contour_inner = new double [len_ctr_inner][2];
    contour_outer = new double [len_ctr_outer][2];


// Add a reference point to close the loop:
    ref_pts_inner[num_ref_pts_inner][0] = ref_pts_inner[0][0];
    ref_pts_inner[num_ref_pts_inner][1] = ref_pts_inner[0][1];
    ref_pts_outer[num_ref_pts_outer][0] = ref_pts_outer[0][0];
    ref_pts_outer[num_ref_pts_outer][1] = ref_pts_outer[0][1];


// Connect Inner Contour
     con_start = 0;
     for(int i=0; i<num_ref_pts_inner; i++)
     { con_end = con_start + resolution_inner[i];
       if(connectors_inner[i]=="line")
       	{ connector_line(i, con_start, con_end, "inner"); }
       
       if(connectors_inner[i]=="circle")
       	{ double xcenter = center_inner[i][0];
          double ycenter = center_inner[i][1];
          connector_circle(xcenter, ycenter, i, con_start, con_end, "inner"); }
       
       if(connectors_inner[i]=="target")
       	{ connector_target_curve(target, targlen, i, con_start, &con_end, "inner"); }
          con_start = con_end; }   
          len_ctr_inner = con_end;
      
       if(num_ref_pts_inner <= 1) len_ctr_inner = 0;


// Connect Outer Contour
     con_start = 0;
     for(int i=0; i<num_ref_pts_outer; i++)
     { con_end = con_start + resolution_outer[i];
       if(connectors_outer[i]=="line")
       	{ connector_line(i, con_start, con_end, "outer"); }
       
       if(connectors_outer[i]=="circle")
       	{ double xcenter = center_outer[i][0];
          double ycenter = center_outer[i][1];
          connector_circle(xcenter, ycenter, i, con_start, con_end, "outer"); }
       
       if(connectors_outer[i]=="target")
       	{ connector_target_curve(target, targlen, i, con_start, &con_end, "outer"); }
        con_start = con_end; }
        len_ctr_outer = con_end;
} // eof


// ====================================================================================
// Get Name
void Part::get_name(string id)
{ partid = id; 
//cout << id << "\n"; 
}

// Get Contour Size
string Part::name() { return partid; }
void Part::get_contour_size(int lenctr, string flag)
	{ if(flag=="inner")
  	{ len_ctr_inner=lenctr; }
  	  else
  	{ len_ctr_outer=lenctr; } }


// ====================================================================================
// Get Contours
void Part::get_contours(double c_in[][2], double c_out[][2],
                        double z_in[],    double z_out[])
{
// ... Inner contour: ...
    contour_inner = new double [len_ctr_inner][2];
    depth_inner   = new double [len_ctr_inner];

  for( int i=0; i<len_ctr_inner; i++)
   { depth_inner[i] = z_in[i];  
  for(int j=0; j<2; j++)
   { contour_inner[i][j] = c_in[i][j]; } }


// ... Outer contour: ...
   contour_outer = new double [len_ctr_outer][2];
   depth_outer   = new double [len_ctr_outer];

  for(int i=0; i<len_ctr_outer; i++)
   { depth_outer[i] = z_out[i];
      
  for(int j=0; j<2; j++)
   { contour_outer[i][j] = c_out[i][j]; } }
}  // eof


// ====================================================================================
// Get Depth
void Part::get_depth(double z_in[], double z_out[])
{

// ... Inner contour: ...
  depth_inner = new double [len_ctr_inner];    
  for( int i=0; i<len_ctr_inner; i++)
   { depth_inner[i] = z_in[i]; }


// ... Outer contour: ...
  depth_outer = new double [len_ctr_outer];
  for( int i=0; i<len_ctr_outer; i++)
   { depth_outer[i] = z_out[i]; }
} // eof


// ====================================================================================
// Spatial Integrals

// Get Volume
// LCC, volumes are calculated for each part according to their contours,
// e.g. "(*BuckCyl).get_volume()".  Is this for a 22.5 deg wedge? (360/16)
// Yes, it only calculates the volume for the unit contour, then geometry revolves it...
void Part::get_volume()
	{ part_volume = 0.0;

// ... integrate over the outer loop ...
    for(int i=0; i<len_ctr_outer-1; i++)
    { double xlo = contour_outer[i][0];
	  double ylo = contour_outer[i][1];
	  double dlo = depth_outer[i];
	  double xhi = contour_outer[i+1][0];
	  double yhi = contour_outer[i+1][1];
	  double dhi = depth_outer[i+1];
      double acf = (yhi-ylo)/(xhi-xlo);
	  double bcf = ylo - acf*xlo;
      double area_loc = (1.0/2.0)*acf*(pow(xhi,2) - pow(xlo,2)) + bcf*(xhi-xlo);
      if(xhi == xlo) area_loc = 0.0;
       part_volume = part_volume + area_loc*(dhi+dlo)/2.0; }


// ... integrate over the inner loop ...
	for(int i=0; i<len_ctr_inner-1; i++)
    { double xlo = contour_inner[i][0];
	  double ylo = contour_inner[i][1];
	  double dlo = depth_inner[i];
	  double xhi = contour_inner[i+1][0];
	  double yhi = contour_inner[i+1][1];
	  double dhi = depth_inner[i+1];
      double acf = (yhi-ylo)/(xhi-xlo);
	  double bcf = ylo - acf*xlo;
      double area_loc = (1.0/2.0)*acf*(pow(xhi,2) - pow(xlo,2)) + bcf*(xhi-xlo);
      if(xhi == xlo) area_loc = 0.0;
       part_volume = part_volume - area_loc*(dhi+dlo)/2.0; }
} // eof

// LCC added to transfer enclosed volumes calculated in G.cpp to DP.cpp
void Part::get_enclosed_volume(double enclosed_volume)
	{ enclvol = enclosed_volume; }

double Part::enclosed_volume()
	{ return enclvol; }
	
	
// ====================================================================================
// Get Volumes
void Part::get_shell_volume(double thickness)
	{ part_volume = surface_area_outer*thickness; }  // eof

double Part::volume()
	{ return part_volume; } // eof

void Part::add_to_volume(double extra) 
	{ part_volume = part_volume + extra; } // eof


// ====================================================================================
// Get Surface Area 
// (LCC not sure about this SA calculation)
// Call by e.g. (*Cryo).surface_area("outer");
double Part::surface_area(string con_id)
	{ if(con_id == "inner") 			return surface_area_inner;
  	  if(con_id == "outer") 			return surface_area_outer;
  	  if(con_id == "plasma_perimeter") 	return perimeter_outer; // for plasma perimeter...
  	  // since function below, line 680 only seems to work for magnets
  	  return 0.0; } // eof


void Part::get_surface_areas() // actual integration
	// ... integrate over the outer loop ...
	{ surface_area_outer = 0.0;
    for(int i=0; i<len_ctr_outer-1; i++) { 
      double xlo = contour_outer[i][0];
	  double ylo = contour_outer[i][1];
	  double dlo = depth_outer[i];
	  double xhi = contour_outer[i+1][0];
	  double yhi = contour_outer[i+1][1];
	  double dhi = depth_outer[i+1];
      double point_distance = sqrt(pow((xhi-xlo),2) + pow((yhi-ylo),2));
      // This "*(dhi+dlo)/2.0" is the midpoint of the line segment
      surface_area_outer    = surface_area_outer + point_distance*(dhi+dlo)/2.0;
      perimeter_outer       = perimeter_outer + point_distance; }
      
// ... integrate over the inner loop ...
    surface_area_inner = 0.0;
    for(int i=0; i<len_ctr_inner-1; i++) { 
      double xlo = contour_inner[i][0];
	  double ylo = contour_inner[i][1];
	  double dlo = depth_inner[i];
	  double xhi = contour_inner[i+1][0];
	  double yhi = contour_inner[i+1][1];
	  double dhi = depth_inner[i+1];
      double point_distance = sqrt(pow((xhi-xlo),2) + pow((yhi-ylo),2));
      surface_area_inner    = surface_area_inner + point_distance*(dhi+dlo)/2.0; }
} // eof


// ====================================================================================
// Cross sections
void Part::get_cross_section(int ncs, double pd, double nd)
{
// This function assumes that the cross section described by an r-z curve is rectangular 
// and perpendicular to the curve. The rectangular cross section is described by: 
   plane_delta[ncs]  = pd; // thickness in r-z plane
   normal_delta[ncs] = nd; // thickness in direction normal to r-z plane
} // eof


// ====================================================================================
// Get Perimeters
// Call e.g. "(*Cryo).perimeter("outer")"
double Part::perimeter(string con_id)
	{ if(con_id == "inner") return perimeter_inner;
 	  if(con_id == "outer") return perimeter_outer;
 	  return 0.0; } // eof

void Part::get_perimeters()
	{ perimeter_outer = 0.0;

// ... integrate over the outer loop ...
    for(int i=0; i<len_ctr_outer-1; i++)
    { double xlo = contour_outer[i][0];
	  double ylo = contour_outer[i][1];
//	  double dlo = depth_outer[i]; // don't know why...b/c it is copied from above and not needed
	  double xhi = contour_outer[i+1][0];
	  double yhi = contour_outer[i+1][1];
//	  double dhi = depth_outer[i+1]; //
      double point_distance = sqrt(pow((xhi-xlo),2) + pow((yhi-ylo),2));
      perimeter_outer       = perimeter_outer + point_distance; }
      perimeter_inner = 0.0;

// ... integrate over the inner loop ...
    for(int i=0; i<len_ctr_inner-1; i++)
    { double xlo = contour_inner[i][0];
	  double ylo = contour_inner[i][1];
//	  double dlo = depth_inner[i]; //
	  double xhi = contour_inner[i+1][0];
	  double yhi = contour_inner[i+1][1];
//	  double dhi = depth_inner[i+1]; //
      double point_distance = sqrt(pow((xhi-xlo),2) + pow((yhi-ylo),2));
      perimeter_inner       = perimeter_inner + point_distance; } 
} // eof


// ====================================================================================
// Get Cross Section Width
double Part::cross_section_width(int ncs)
	{ return plane_delta[ncs]; }  // thickness in r-z plane

// Get Cross Section Depth
double Part::cross_section_depth(int ncs)
	{ return normal_delta[ncs]; }  // thickness normal to r-z plane (poloidal direction)



// ====================================================================================
// ====================================================================================
// Material and Cost Operators

void Part::get_number_of_material_components(int num)
{ len_material = num; } // eof


// Setup subroutine, these are called from Cost.cpp to cost each component
// (1) material name, (2) volume fraction, (3) density fraction, and (4) the cost base.
void Part::get_material_component(string matname, double vfrac, double denfrac, int matn)						  
{    material_name[matn]           = matname;
	 // add description, purpose
     volume_fraction[matn]         = vfrac;
     density_fraction[matn]        = denfrac;
     material_density[matn]        = density(matname);
     material_cost_unit_mass[matn] = cost_unit_mass(matname); // RHS is a function
     material_cost_base[matn]      = cost_base(matname); 
//     material_description[matn]    = description(matname); // LCC added
//cout << "matname = " << matname << "\n";
//cout << "description = " << description(matname) << "\n";

// Printout
	string filedir, filepath, extension; filedir = "./OutputParts/"; extension = ".data";
    filepath = filedir + "zPartCompositions" + extension; ofstream fileout(filepath.c_str(), ios::app);
// TO DO: if file already exists, delete it. Also add header

// Print all material compositions in "zPartCompositions" as given in Cost.cpp
fileout << name() << "	;";
fileout	<< matname << "	;"
//		<< description(matname) << "	;"
 		<< vfrac << "	;"
		<< denfrac << "	;"
 		<< density(matname) << "	;"
 		<< cost_unit_mass(matname) << "	;"
 		<< cost_base(matname) << "\n";
} // eof


// ****TESTING - trying to get material descriptions and purposes from input file to print out
/*string Part::description(string matname)
	{ string desc = " ";
	// Read material properties from materials.data input file
	ifstream fin_mat;  fin_mat.open("InputData/materials.data");
	fin_mat.width(10); fin_mat.precision(2);
	fin_mat.setf(ios_base::uppercase, ios_base::scientific);
	int next = 60;   // next line of materials.data after * at end
	string description[22];

for(int i=0, j=0; i<22; i++)
    { std::ostringstream stm;
      fin_mat.seekg((j+1)*next); fin_mat >> description[i]; // description
      stm << description[i];
      desc += stm.str() + " " + material_name[i] + " ";
      j=j+7;
      cout << "description = " << desc << "\n";	  
      return desc; // function returns this compiled string
      }
} */

/*
// Get material description and purpose
string Part::description(string matname)						  
{	string desc[22] = " ";
	// Read material properties from materials.data input file
	ifstream fin_mat;  fin_mat.open("InputData/materials.data"); // 1.5 of 3
	fin_mat.width(10); fin_mat.precision(2);
	fin_mat.setf(ios_base::uppercase, ios_base::scientific);
	int next = 60;   // next line of materials.data after * at end
	string name[22]; // should make this a global variable in case more materials are added
	string description[22];
	string purpose[22];

for(int i=0, j=0; i<22; i++)
    { fin_mat.seekg(j*next);     fin_mat >> name[i];        // name
	  fin_mat.seekg((j+1)*next); fin_mat >> description[i]; // description
	  fin_mat.seekg((j+2)*next); fin_mat >> purpose[i];     // purpose
	  //description[i] = prop[c+2]; // put description in a separate vector to use below
	  //purpose[i]     = prop[c+3]; // put purpose in a separate vector to use below
      j=j+7;
      } // iterate by 5 material properties & 7 lines per material name
//cout << "description = " << description[1] << "\n";
//cout << "purpose     = " << purpose[1] << "\n";             

// Find material description and purpose when name matches and return it
	for(int i=0; i<22; i++)
    	{ if(matname == name[i]) return description[i]; 
cout << "description = " << description[i] << "\n";
    	}

//return desc;
fin_mat.close();
} // eof
*/

// ====================================================================================
// Density of Materials [kg/m3]

double Part::density(string matname) {
// *** Print-out for displaying material properties was moved to Aries.cpp
// This is not the correct way to get these properties since this function is called
// multiple times...

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

for(int i=0, j=0, c=0; i<22; i++) // loop through all 22 materials
    { fin_mat.seekg(j*next);     fin_mat >> name[i];        // name
//      fin_mat.seekg((j+1)*next); fin_mat >> description[i]; // description
//      fin_mat.seekg((j+2)*next); fin_mat >> purpose[i];     // purpose
      fin_mat.seekg((j+3)*next); fin_mat >> prop[c];        // density
	  density[i] = prop[c]; // put density in a separate vector to use below
//cout << "prop[i] = " << density[i] << "\n";
      c=c+5; j=j+7;
      } // iterate by 5 material properties & 7 lines per material name

// Find density when name matches and return it
	for(int i=0; i<22; i++)
    { if(matname == name[i]) return density[i]; }

return 0.0; // if nothing in density[i] is returned
fin_mat.close();
} // eof


// ====================================================================================
// Liquid Metal Mass [kg]
// written by ZD after email conversation with LW ~1/30/10
// this is for 1/16th of part since part_volume is 1/16th
double Part::liquid_metal_mass()
	{ mass_of_liquid_metal = 0.0;
   	for(int i=0; i<50; i++)
   		{ if(material_name[i] == "PbLi")
        { mass_of_liquid_metal = part_volume*material_density[i]*volume_fraction[i]; 
//    cout << "LM Mass = " << mass_of_liquid_metal << " = " << part_volume << " * " 
//    					 << material_density[i] << " * " << volume_fraction[i] << "\n";
    } }   
return mass_of_liquid_metal; 
} // eof


// ====================================================================================
// Unit Cost Mass [$/kg]

double Part::cost_unit_mass(string matname) 
{
// Read material properties from materials.data input file
	ifstream fin_mat;  fin_mat.open("InputData/materials.data"); // 2 of 3
	fin_mat.width(10); fin_mat.precision(2);
	fin_mat.setf(ios_base::uppercase, ios_base::scientific);
	int next = 60;  // next line
	string name[22]; double prop[131]; double unitcost[22];

for(int i=0, j=0, c=0; i<22; i++)
    { fin_mat.seekg(j*next);     fin_mat >> name[i];   // name
	  fin_mat.seekg((j+4)*next); fin_mat >> prop[c+1]; // cost/kg
	  unitcost[i] = prop[c+1]; // put unit cost in a separate vector to use below
      c=c+5; j=j+7; } // iterate by 5 material properties per name
             
// Find unit cost when name matches and return it
	for(int i=0; i<22; i++)
    	{ if(matname == name[i]) return unitcost[i]; }
      
return 0.0; // if nothing is returned by unitcost
fin_mat.close();
} // eof


// ====================================================================================
// Material Properties Table

// This is now in input file "materials.data" and has updated descriptions and properties
// Data is from 'Fabrication Costs per kg Rev# (date).xls' from L. Waganer
// Also, magnet properties taken from ASC Docs "PF,TFMagnet.doc" from L. Bromberg

// Name				"name" in LW spreadsheet 		density kg/m3, cost_unit_mass $/kg, cost base year
//--------------------------------------------------------------------------------------
// 1  B-FS     		= "FS, RA, Blanket structure" borated			7800	35.54	2009 for blanket composition, was 37.42			 	        
// 2  PbLi     		= "Li17Pb83, liquid (Enriched 90% Li6) @650C"	8897	8.79	2009 must input in CA.cpp, ~1238 in 27.1	
// 3  RAFS      	= "FS, RA, Large, Complex Struct (composite)"	7800	43.44	2009 for VV, use composite to include eveything. Laila - RAFS is generic name for F82H (Japanese) and Eurofer (European steel), she suggests 2x for safety compliance	
// 4  SiC-comp 		= "SiC/SiC composite structure"					3200	584.64	2009		
// 5  SS       		= "304 SS very large, simple, accurate"			7800	81.46	2009 use for support structure
// 6  Void     		= "void"										0.0		0.0		2009	
// 7  W        		= "W, cast"										19300	120.37	2009
// 8  WC       		= "WC, balls, plates"							15600	34.39	2009
// 9  TF-coil      	= NEED composition								7500	50.912	1992
// 10 PF-coil  		= NEED composition								7500	82.919	1992
// 11 Nb3Sn       	= PFMagnet - "LT Supercond"						7500	92.70	1992
// 12 Cu       		= TFMagnet - "Copper"							7300	65.4	1992
// 13 YBCO     		= TFMagnet - "HT Supercond"						6200	92.7	1992
// 14 Ag       		= NEED, data from where?						10500	201.1	1992
// 15 CeO2       	= NEED, data from where?						7130	16.516	1992
// 16 Polyamide    	= TFMagnet - "Polyamide"						1920	16.516	1992
// 17 Inconel		= TFMagnet - "Inconel"							8440	96.02	2009 was 50.912 1992, changed to LW's table
// 18 JK2LB       	= "JK2LB steel (Conv. Fab, Tokamak Coils)"		8000	95.0	2009 
// 19 SS316			= TFMagnet - Stainless Steel 316				7860	68.657	1992 was "structconserv" 8000, 100.0, 1992	
// 20 FS      		= "FS, RA, Blanket structure,"					7800	89.42	2009
// 21 He       		= "99.999% pure, Qingdao Baigong.. China (LW)	0.1786	70.0*	2011
// 22 RAFS-div		= "FS, RA, divertor structure (25% surcharge)"	7800	111.77	2009 for div
// 22 RAFS-cryost	= "FS, RA, Mod size, Welded Dbl-Walled Assy"	7800	28.06	2009 for cryostat
//--------------------------------------------------------------------------------------
// * Helium is zero in the table since it is added at startup and not part of the building contract
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



// ====================================================================================
// Cost Base [year]

int Part::cost_base(string matname) {

// Read material properties from materials.data input file
	ifstream fin_mat;  fin_mat.open("InputData/materials.data"); // 3 of 3
	fin_mat.width(10); fin_mat.precision(2);
	fin_mat.setf(ios_base::uppercase, ios_base::scientific);
	int next = 60;  // next line
	string name[22]; double prop[131]; int costbase[22];

for(int i=0, j=0, c=0; i<22; i++)
    { fin_mat.seekg(j*next);     fin_mat >> name[i];   // name
	  fin_mat.seekg((j+5)*next); fin_mat >> prop[c+2]; // cost base year
	  costbase[i] = prop[c+2]; // put cost base in a separate vector to use below
      c=c+5; j=j+7; } // iterate by 5 material properties per name

// Find cost base year when name matches and return it
	for(int i=0; i<22; i++)
    	{ if(matname == name[i]) return costbase[i]; } 

return 0; // if nothing is returned by costbase
fin_mat.close();
} // eof


// ====================================================================================
// Cost of Part [$]
// called from Cost.cpp once specific material properties are taken from tables
void Part::get_cost(int baseYear)
{  total_cost = 0.0;
   for(int i=0; i<len_material; i++)
   { int locDefYear = material_cost_base[i];
     total_cost = total_cost +
   				  part_volume*
   				  volume_fraction[i]*
   				  density_fraction[i]*
   				  material_density[i]*
   				  material_cost_unit_mass[i]*
   				  inflation(locDefYear, baseYear); // returns inflation factor
//cout << "locDefYear = " << locDefYear << "\n";
//cout << "baseYear = " << baseYear << "\n";

} } // eof


// ====================================================================================
// Mass of Parts [kg]
void Part::get_mass()
{	total_mass = 0.0;
   	for(int i=0; i<len_material; i++) // len_material = # material components
   	{ total_mass = total_mass +
                   part_volume*
                   volume_fraction[i]*
                   density_fraction[i]*
                   material_density[i]; } } // eof


// ====================================================================================
// Sub Costs
void Part::get_sub_cost(int count, double cost) { partial_cost[count] = cost; } // eof

// Functions that return values to other modules
void   Part::add_to_cost(double extra) { total_cost = total_cost + extra; } // eof
void   Part::add_to_mass(double extra) { total_mass = total_mass + extra; } // eof
double Part::cost() 				   { return total_cost; } 				// eof
double Part::sub_cost(int count)	   { return partial_cost[count]; } 		// eof
double Part::mass() 				   { return total_mass; } 				// eof
int    Part::number_material_components() { return len_material; } 		    // eof
// also (*Part).cost_base("He")


// ====================================================================================
// Part Material Composition
string Part::material_composition()
	{ string output = " ";
    for(int i=0; i<len_material; i++)
   	{ std::ostringstream stm;
      stm << volume_fraction[i]*100.0;
      output += stm.str() + "% " + material_name[i] + " "; } // += adds to string
	  return output; // function returns this compiled string
}


// List Part Components
string Part::part_components(int num)
// num is part material #, use (*Part).number_material_components() to get # components
	{ string output = " ";
   	  std::ostringstream stma, stmb, stmc, stmd, stme;
      stma << volume_fraction[num]; // required for putting numbers into a string
      stmb << density_fraction[num];
      stmc << material_density[num];
      stmd << material_cost_unit_mass[num];
      string matname = material_name[num];
      stme << cost_base(matname); // cost_base input must be the material name string

// material_name is a string, density & volume_fraction are doubles
      output = matname 		+ "	;" + 
      		   stma.str() 	+ "	;" + 
      		   stmb.str() 	+ "	;" + 
      		   stmc.str() 	+ "	;" +
      		   stmd.str() 	+ "	;" +
      		   stme.str()	;
//cout << "desc = " << description[num] << "\n";
	  return output; // function returns this single string
}


// Full Part Printout
string Part::part_printout()
	{ string output = " ";
	for(int i=0; i<len_material; i++)
	{ std::ostringstream stma, stmb, stmc, stmd, stme, stmf, stmg;
	  stma << part_components(i);
	  if(i == 0) {
	  stmc	<< part_volume*16.0;
 	  stmd	<< surface_area("outer")*16.0/2.0;
	  stme	<< mass_of_liquid_metal*16.0/1000.0;
 	  stmf	<< total_mass;
 	  stmg	<< total_cost*16.0/1.e6;
 	  output += "	;" 		+
 	  			stma.str() 	+ "	;" + 
// 	  			stmb.str() 	+ "	;" +
 	  			stmc.str() 	+ "	;" +
				stmd.str() 	+ "	;" +
      		    stme.str() 	+ "	;" +
      		    stmf.str() 	+ "	;" +
      		    stmg.str()	+ "\n" ; }	  
 	  else { output += "	;" + stma.str() + "	\n"; }
	}
return output;}

// ====================================================================================
// ====================================================================================
// Inflation Factors 

// Also in CostingAccount.cpp line ~1945 at bottom and prints out in Aries.cpp
double Part::inflation(int yearFrom, int yearTo)  // Year input from input file costing.data
{
// NEW - Read GDP inflation factors from inflationfactors.data input file
	ifstream fin_mat;  fin_mat.open("InputData/inflationfactors.data");
	fin_mat.width(10); fin_mat.precision(2);
	fin_mat.setf(ios_base::uppercase, ios_base::scientific);
	int next = 15;  // next line
	int year[42]; double factor[85];

// Read the year data from the input file
for(int i=3, j=0; i<86; j++)
    { fin_mat.seekg(i*next);     fin_mat >> year[j];   // year
      fin_mat.seekg((i+1)*next); fin_mat >> factor[j]; // factor
      i=i+2; 
//cout << year[j] << " " << factor[j] << "\n";
}

// Calculate ratio from given input years
	double cf_from, cf_to;
	for(int i=0; i<42; i++)
		{ if(yearFrom == year[i]) cf_from = factor[i]; 
	 	  if(yearTo   == year[i]) cf_to   = factor[i]; }
//cout << "ratio = " << cf_from/cf_to << "\n";

return cf_to/cf_from; // inflation() function returns this ratio
fin_mat.close();
} //eof



// ====================================================================================
// ====================================================================================
// Display Data, print IBL files of part definitions for CAD drawings

void Part::show() // called by DP ~ 1060
{  string filedir, filepath, filepathIBL, extension, extensionIBL;
   filedir 		= "./OutputParts/"; extension = ".data"; extensionIBL = ".ibl";
   filepath 	= filedir + partid + extension;
   filepathIBL 	= filedir + partid + extensionIBL;
   ofstream fileout(filepath.c_str()); ofstream fileoutIBL(filepathIBL.c_str());
   double len = (len_ctr_inner > len_ctr_outer) ? len_ctr_inner : len_ctr_outer; 

// Not needed:
//cout << "\n" << "Data for part " << partid << " available in " << filepath << "\n";
//cout << partid << "\n";

// File *.data for processing by Matlab
    for(int i=0; i<len; i++)
    { if(i < len_ctr_inner)
    { fileout << contour_inner[i][0] << "	"					<< contour_inner[i][1] << "	"; }
    else
    { fileout << contour_inner[len_ctr_inner-1][0] << "	" 		<< contour_inner[len_ctr_inner-1][1] << "	"; }
    
    if(i < len_ctr_outer)
    { fileout << contour_outer[i][0] << "	" 					<< contour_outer[i][1] << "	"; }
    else
    { fileout << contour_outer[len_ctr_outer-1][0] << "	" 		<< contour_outer[len_ctr_outer-1][1] << "	"; }
    
    if(i < num_ref_pts_inner)
    { fileout << ref_pts_inner[i][0] << "	" 					<< ref_pts_inner[i][1] << "	"; }
    else
    { fileout << ref_pts_inner[num_ref_pts_inner-1][0] << "	" 	<< ref_pts_inner[num_ref_pts_inner-1][1] << "	"; }

    if(i < num_ref_pts_outer)
    { fileout << ref_pts_outer[i][0] << "	" 					<< ref_pts_outer[i][1]; }
    else
    { fileout << ref_pts_outer[num_ref_pts_outer-1][0] << "	" 	<< ref_pts_outer[num_ref_pts_outer-1][1]; }
    fileout << "\n"; }


// File *.ibl for processing by PRO/Engineer
const double pi = 4.0*atan(1.0); // pi = 3.14159....
int mpcs = 72; // number of plasma cross sections, created by rotating cross-section every 4.5 deg

	if(partid != "Plasma_DCLL") 
	{ fileoutIBL << "/* Number of cross-section = 1 */\n"; }
	else
	{ fileoutIBL << "/* Number of period = 2 */\n"; }

	if(partid == "Separatrix") 
	{ fileoutIBL << "open\n";
	  fileoutIBL << "arclength\n\n"; } 
	else
	{ fileoutIBL << "closed\n";
	  fileoutIBL << "pointwise\n\n"; }
	  
	  
// Treatment of any entity other than "plasma"
	if(partid != "Plasma_DCLL") 
	{ fileoutIBL << "begin section ! 1\n";
	  fileoutIBL << "begin curve ! 1\n";
	  int kpoint; // data point number in the output file


// Write the coordinates of the outer curve:
kpoint = 0;
for(int i=0; i<len_ctr_outer; i++)
{	kpoint++;
	double xout = contour_outer[i][0]*cos(0.0);
	double yout = contour_outer[i][0]*sin(0.0);
	double zout = contour_outer[i][1];

 	fileoutIBL << setw(12) << kpoint;
 	fileoutIBL << setw(15) << setprecision(5) << xout;
 	fileoutIBL << setw(15) << setprecision(5) << yout;
 	fileoutIBL << setw(15) << setprecision(5) << zout << "\n";

// TF Magnet
if(partid == "TFMagnet" & (i==0 | i==len_ctr_outer-2) )
{  double xoutplus = contour_outer[i+1][0]*cos(0.0);
   double youtplus = contour_outer[i+1][0]*sin(0.0);
   double zoutplus = contour_outer[i+1][1];
   
 // Add 9 data points in the flat portion of the TF magnet for Pro/E processing
for(int j=0; j<9; j++)
{  kpoint++;
   double xoutins = xout + (j+1)*(xoutplus-xout)/10.0;
   double youtins = yout + (j+1)*(youtplus-yout)/10.0;
   double zoutins = zout + (j+1)*(zoutplus-zout)/10.0;

   fileoutIBL << setw(12) << kpoint;
   fileoutIBL << setw(15) << setprecision(5) << xoutins;
   fileoutIBL << setw(15) << setprecision(5) << youtins;
   fileoutIBL << setw(15) << setprecision(5) << zoutins << "\n"; } } }


// Separatrix
if(partid != "Separatrix")
{	double xout = contour_outer[0][0]*cos(0.0);
	double yout = contour_outer[0][0]*sin(0.0);
	double zout = contour_outer[0][1];

	fileoutIBL << setw(12) << kpoint + 1;
	fileoutIBL << setw(15) << setprecision(5) << xout; 
	fileoutIBL << setw(15) << setprecision(5) << yout;
	fileoutIBL << setw(15) << setprecision(5) << zout << "\n\n"; }
else
	{ fileoutIBL << "\n"; }
	
//	fileoutIBL << "begin section ! 2\n";
//	fileoutIBL << "begin curve ! 2\n";


// Write the coordinates of the inner curve:
kpoint = 0;
for(int i=0; i<len_ctr_inner; i = i++)
{ 	kpoint++;
 	double xout = contour_inner[i][0]*cos(0.0);
 	double yout = contour_inner[i][0]*sin(0.0);
 	double zout = contour_inner[i][1];

 	fileoutIBL << setw(12) << kpoint;
 	fileoutIBL << setw(15) << setprecision(5) << xout;
 	fileoutIBL << setw(15) << setprecision(5) << yout;
	fileoutIBL << setw(15) << setprecision(5) << zout << "\n";


// TF Magnet (replicated?)
if(partid == "TFMagnet" & (i==0 | i==len_ctr_outer-2) )
{  double xoutplus = contour_inner[i+1][0]*cos(0.0);
   double youtplus = contour_inner[i+1][0]*sin(0.0);
   double zoutplus = contour_inner[i+1][1];

// add 9 data points in the flat portion of the TF magnet for Pro/E processing
for(int j=0; j<9; j++) 
{  kpoint++;
   double xoutins = xout + (j+1)*(xoutplus-xout)/10.0;
   double youtins = yout + (j+1)*(youtplus-yout)/10.0;
   double zoutins = zout + (j+1)*(zoutplus-zout)/10.0;
//cout << xoutins << "\n";
   
   fileoutIBL << setw(12) << kpoint;
   fileoutIBL << setw(15) << setprecision(5) << xoutins;
   fileoutIBL << setw(15) << setprecision(5) << youtins;
   fileoutIBL << setw(15) << setprecision(5) << zoutins << "\n"; } } }


// Separatrix (replicated?)
if(partid != "Separatrix")
{	double xout = contour_inner[0][0]*cos(0.0);
	double yout = contour_inner[0][0]*sin(0.0);
	double zout = contour_inner[0][1];

	fileoutIBL << setw(12) << kpoint + 1;
	fileoutIBL << setw(15) << setprecision(5) << xout; 
	fileoutIBL << setw(15) << setprecision(5) << yout;
	fileoutIBL << setw(15) << setprecision(5) << zout << "\n\n"; }
else
	{ fileoutIBL << "\n"; } }


else // anything other than Plasma_DCLL

// Special treatment for plasma entity
{ int npts_max = 101; // 101 = maximum number of data points allowed in the output file.
                      // The code will skip len_ctr_outer/npts_max + 1 data points
                      // between each two points printed in the output file.
                      // This is done in order to make the IBL output file smaller.

for(int j=0; j < mpcs; j++) // loop over cross sections sliced by 5 degrees in theta direction
{ fileoutIBL << " begin section ! " << j+1 << "\n";
  fileoutIBL << " begin curve ! 1\n";


// Write the coordinates of the outer curve:
int kpoint = 0; // data point number in the output file
for(int i=0; i<len_ctr_outer; i = i + len_ctr_outer/npts_max +1)
{	kpoint++;
	double xout = contour_outer[i][0]*cos(j*2.0*pi/mpcs);
 	double yout = contour_outer[i][0]*sin(j*2.0*pi/mpcs);
	double zout = contour_outer[i][1];

 	fileoutIBL << setw(12) << kpoint;
	fileoutIBL << setw(15) << setprecision(5) << xout;
 	fileoutIBL << setw(15) << setprecision(5) << yout;
 	fileoutIBL << setw(15) << setprecision(5) << zout << "\n"; }

 	double xout = contour_outer[0][0]*cos(j*2.0*pi/mpcs);
 	double yout = contour_outer[0][0]*sin(j*2.0*pi/mpcs);
 	double zout = contour_outer[0][1];

 	fileoutIBL << setw(12) << kpoint + 1;
 	fileoutIBL << setw(15) << setprecision(5) << xout; 
 	fileoutIBL << setw(15) << setprecision(5) << yout;
 	fileoutIBL << setw(15) << setprecision(5) << zout << "\n\n";

} // end loop over cross sections
} // end special treatment for plasma
  // fileoutIBL << "/* end of file */";

}  // eof