# reads in ribs2gpdata.sh from delme.txt , then using ~greg/smooth.r to interpolate the ribs
# , upping the res of the data points.

# does the same for the keel line, to include the bow shape.

# then loops from rib to rib, picking making surfaces out of the same z and the next one down (or up??)
# starts at rib 1 and looks towards the front, so that the bow can be interpreted as rib 0.

# To me they are ribs, the boat plan called them stations.

# data is in boat plan inches, 302 " long, 72 " wide.

# R code cannot be told how many points to create, but bigger input numbers leads to more output
# interpolitad points, so convert to full size boat mm.

# As the R code doesn't support 3d coords, had to write my own. Alas it has bugs, as the maths
# doesnt really support s[-1] through to s[+1]. 0 to 1 was supported, but -1 support broke anything not -1.

# bash ribs2gpdata.sh ; bash gpdata2obj.sh 

gawk -- '{
   # reset state with any major header, upper/lower/keel
   if( $0 ~ /^# / )
      s="";

   # set rib state with any minor header, ie, #rib <rib number>
   if( (s=="upper" || s=="lower") && match($0,"^#rib") ){ n=0; rib=$2; num_ribs=rib; printf "set rib to %s\n",rib; }

   #printf "state=%s, rib=%s, n=%s\n",s,rib,n;

   # store the numbers
   scale=25.4*1;  #(170/1800);  # inches to mm , and then final scaling, making the maximum beam (width) 170 mm, small enough to fit on the build plate.
   scale_abs=1;  #170/1800;  # and repeat for the all magic numbers...
   # ***except*** that uniform scaling had side effects, eaiser to scale downstream, there are only boat.obj and block.obj coming from this code.

   morescale=1;  # make the magnitude even bigger, to get more data from interpolation. Doesnt do anything, same number of output points
   if( match($1,"^[0-9-]")>0 )
    {
      if( s=="upper" )
       { upp_x[rib,n]=$1*scale; upp_y[rib,n]=$2*scale; upp_z[rib,n++]=$3*scale;  upp_num[rib]=n; }
      else if( s=="lower" )
       { low_x[rib,n]=$1*scale; low_y[rib,n]=$2*scale; low_z[rib,n++]=$3*scale;  low_num[rib]=n; }
      else if( s=="keel" )
       { keel_x[n]=$1*scale; keel_y[n]=$2*scale; keel_z[n++]=$3*scale;  keel_num=n; }
      else if( s=="shaft" )
       { shaft_x[n]=$1*scale; shaft_y[n]=$2*scale; shaft_z[n++]=$3*scale;  shaft_num=n; }
      else if( s=="block" )
       { block_x[n]=$1*scale; block_y[n]=$2*scale; block_z[n++]=$3*scale;  block_num=n; }
    }


   # set major state as required, upper, lower ribs, or keel.
   if( $0 == "# upper ribs" )
    { s="upper"; n=0; rib=0; }
   else if( $0 == "# lower ribs" )
    { s="lower"; n=0; rib=0; }
   else if( $0 == "# keel" )
    { s="keel"; n=0; }
   else if( $0 == "# shaft" )
    { s="shaft"; n=0; }
   else if( $0 == "# block" )
    { s="block"; n=0; }

 } function inter( y1, z1,  y2, z2,  zt )
 {
   return (y2-y1)/(z2-z1) * (zt-z1) + y1;
 } END {

   ### prop supporting block, as a sep object to import into openscad
   if(1)
    {
      fn="block.obj";
      printf "# block supporing prop shaft\n" > fn;
      vertex=1;
      topfront_x=block_x[3];   topfront_y=0;  topfront_z=block_z[3];
      topback_x=block_x[0];    topback_y=0;   topback_z=block_z[0];
      botfront_x=block_x[2];   botfront_y=0;  botfront_z=block_z[2];
      botback_x=block_x[1];    botback_y=0;   botback_z=block_z[1];
      d2r=atan2(1,1)/45;
      da=2;
      for(a=0;a<360;a+=da)
       {
         centertop_x=(topfront_x+topback_x)/2;   disttop_w=topfront_x-centertop_x;  disttop_y=disttop_w*0.2;
         centerbot_x=(botfront_x+botback_x)/2;   distbot_w=botfront_x-centerbot_x;  distbot_y=distbot_w*0.2;

         bt_x[a]=centertop_x + cos(a*d2r)*disttop_w;  bt_y[a]=sin(a*d2r)*disttop_y;
         printf "v %f %f %f\n",bt_x[a], bt_y[a], topfront_z >> fn;
         bt_v[a]=vertex++;

         bb_x[a]=centerbot_x + cos(a*d2r)*distbot_w;  bb_y[a]=sin(a*d2r)*distbot_y;
         printf "v %f %f %f\n",bb_x[a], bb_y[a], botfront_z >> fn;
         bb_v[a]=vertex++;
       }
      printf "v %f %f %f\n", centertop_x, 0, topfront_z >> fn;
      bt_center_v=vertex++;
      printf "v %f %f %f\n", centerbot_x, 0, botfront_z >> fn;
      bb_center_v=vertex++;

      for(a=0;a<360;a+=da)
       {
         dn=a+da;  if(dn>=360)dn-=360;
         printf "f %d %d %d\n",bt_v[a], bb_v[a], bt_v[dn] >> fn;   # sides
         printf "f %d %d %d\n",bb_v[a], bb_v[dn], bt_v[dn] >> fn;

         printf "f %d %d %d\n",bt_center_v, bt_v[a], bt_v[dn] >> fn;   # top (wide) end
         printf "f %d %d %d\n",bb_center_v, bb_v[dn], bb_v[a] >> fn;   # bottom (narrower) end
       }

      close(fn);
    }

   #for(i=0;i<block_num;i++)
   # {
   #   printf "%d\t%d %d %d\n",i,block_x[i], block_y[i], block_z[i];
   # }

   # shaft maths
   for(i=0;i<shaft_num;i++)
      printf "%d\t%d %d %d\n",i,shaft_x[i], shaft_y[i], shaft_z[i];

   printf "shaft angle of %f from hozi\n",180 - atan2( shaft_z[0]-shaft_z[1], shaft_x[0]-shaft_x[1] )/(atan2(1,1)/45);
   printf "shaft length of %f\n",sqrt( (shaft_z[0]-shaft_z[1])^2 + (shaft_x[0]-shaft_x[1])^2 );

   # shaft coards are 3068,341 (engine end) to (7675,-89) (just past prop nut).
   # => angle is.. 5.340021 degrees from horiz and 4626.94 mm long.
   # where the length should then be scaled by 170/1800, as should boat.obj and block.obj

   ### main hull code

   for(r=1;r<=num_ribs;r++)
    {
      printf "" > "test.dat";
      for(p=0;p<upp_num[r];p++)
         printf "%f\t%f\t-1\n",upp_y[r,p],upp_z[r,p] >> "test.dat";
      close("test.dat");
      system("R --no-save < ~/smooth.r > /dev/null 2> /dev/null");

      p=0;
      x=upp_x[r,1];
      while( (getline < "delme.csv")==1 )
       {
         upp_x[r,p]=x;  upp_y[r,p]=$1;  upp_z[r,p++]=$2;
         upp_num[r]=p;         
         #printf "upp %d,%d  %f,%f,%f\n",r,p, upp_x[r,p-1], upp_y[r,p-1], upp_z[r,p-1];
       }
      close("delme.csv");

    }


   for(r=1;r<=num_ribs;r++)
    {
      printf "" > "test.dat";
      for(p=0;p<low_num[r];p++)
         printf "%f\t%f\t-1\n",low_y[r,p],low_z[r,p] >> "test.dat";
      close("test.dat");
      system("R --no-save < ~/smooth.r > /dev/null 2> /dev/null");

      p=0;
      x=low_x[r,1];
      while( (getline < "delme.csv")==1 )
       {
         low_x[r,p]=x;  low_y[r,p]=$1;  low_z[r,p++]=$2;
         low_num[r]=p;         

         #printf "low %d,%d  %f,%f,%f\n",r,p, low_x[r,p-1], low_y[r,p-1], low_z[r,p-1];
       }
      close("delme.csv");
    }

   printf "" > "test.dat";
   for(p=0;p<keel_num;p++)
      printf "%f\t%f\t-1\n",keel_x[p],keel_z[p] >> "test.dat";
   close("test.dat");
   system("R --no-save < ~/smooth.r > /dev/null 2> /dev/null");

   p=0;
   y=keel_y[1];
   while( (getline < "delme.csv")==1 )
    {
      keel_x[p]=$1;  keel_y[p]=y;  keel_z[p++]=$2;
      keel_num=p;         
    }
   close("delme.csv");

   # have now got the upper and lower ribs, and the keel section in a higher res.
   # before tracing, I feel itll be better to rotate around x, so that upper and lower ribs can be combined.


   # merge upper and lower rib coords
   # dont merge, it would now be better seperate
   if(0)
   for(r=1;r<=num_ribs;r++)
    {
      p2=upp_num[r];
      for(p=0;p<low_num[r];p++)
       {
         upp_x[r,p2]=low_x[r,p];  upp_y[r,p2]=low_y[r,p];  upp_z[r,p2++]=low_z[r,p];
       }
      upp_num[r]=p2;
    }

   # rotate all coords by 45 degrees
   # dont rotate, now we can do 3D splines, there is no need. And besides, rotating didnt help.
   if(0)
    {
      a=atan2(1,1);
      for(p in keel_z)
       {
         d=sqrt( keel_y[p]**2 + keel_z[p]**2 );
         y=keel_y[p];
         keel_y[p]=cos(a)*keel_y[p] - sin(a)*keel_z[p];
         keel_z[p]=sin(a)*y         + cos(a)*keel_z[p];  
         #printf "dist=%f, %f\n", sqrt( keel_y[p]**2 + keel_z[p]**2 ) - d, d;
       }
      for(r=1;r<=num_ribs;r++)
         for(p=0;p<upp_num[r];p++)
          {
            d=sqrt( upp_y[r,p]**2 + upp_z[r,p]**2 );
            y=upp_y[r,p];
            upp_y[r,p]=cos(a)*upp_y[r,p] - sin(a)*upp_z[r,p];
            upp_z[r,p]=sin(a)*y         + cos(a)*upp_z[r,p];
            #printf "%d,%d dist=%f, %f\n", r,p, sqrt( upp_y[r,p]**2 + upp_z[r,p]**2 ) - d, d;
          }
    }

   # extract cooards of constant Z (height), from bow to stern, or at least as far as poss - bow is higher than stern, together are
   # concave. Now they are rotated by 45 degrees, hard to say what shape it is...


   # obj header and ribs and keel . Just FYI, faces dont refer to these ribs or keel (yet, will need join points)
   fn="boat.obj";
   printf "# gpdata2obj.sh\n#\n# object Boat Hull\n#\n" > fn;
   # save rib vertices to obj, to see where the control points are in the .xfig files.
   vertex=1;  # track vertex number, so the faces know which vertices to refer to.
   for(r=1;r<=num_ribs;r++)
    {
      for(p=0;p<upp_num[r];p++)
       {
         printf "v %f %f %f\n",upp_x[r,p],upp_y[r,p],upp_z[r,p] >> fn;  # original
         upp_v[r,p]=vertex++;
         printf "v %f %f %f\n",upp_x[r,p],-upp_y[r,p],upp_z[r,p] >> fn; # mirror
         vertex++;
       }

      for(p=0;p<low_num[r];p++)
       {
         printf "v %f %f %f\n",low_x[r,p],low_y[r,p],low_z[r,p] >> fn;  # original
         low_v[r,p]=vertex++;
         printf "v %f %f %f\n",low_x[r,p],-low_y[r,p],low_z[r,p] >> fn;  # mirror
         vertex++;
       }
    }
   for(p in keel_z)
    {
      printf "v %f %f %f\n", keel_x[p], keel_y[p], keel_z[p] >> fn;
      keel_v[p]=vertex++;
    }

   #num_ribs=3;  # force rib count, to concentrate on just the bow

   # save out ribs and keel to .csv file
   printf "" > "delme2.txt";
   for(p in keel_z)
      printf "%f %f %f\n",keel_x[p],keel_y[p],keel_z[p] >> "delme2.txt";
   for(r=1;r<=num_ribs;r++)
      for(p=0;p<upp_num[r];p++)
         printf "%f %f %f\n",upp_x[r,p],upp_y[r,p],upp_z[r,p] >> "delme2.txt";
   for(r=1;r<=num_ribs;r++)
      for(p=0;p<low_num[r];p++)
         printf "%f %f %f\n",low_x[r,p],low_y[r,p],low_z[r,p] >> "delme2.txt";

   # starting with top of upper rib, move to bottom of upper rib by this fraction
   # so need to calculate total rib length. Do same for lower rib too.
   for(r=1;r<=num_ribs;r++)
    {
      upp_len[r,0]=0;
      lastx=upp_x[r,0];  lasty=upp_y[r,0];  lastz=upp_z[r,0];
      for(p=1;p<upp_num[r];p++)
       {
         upp_len[r,p]=upp_len[r,p-1] + sqrt( (lastx-upp_x[r,p])**2 + (lasty-upp_y[r,p])**2 + (lastz-upp_z[r,p])**2 );
         lastx=upp_x[r,p];  lasty=upp_y[r,p];  lastz=upp_z[r,p];
       }
      upp_totlen[r]=upp_len[r,p-1];

      low_len[r,0]=0;
      lastx=low_x[r,0];  lasty=low_y[r,0];  lastz=low_z[r,0];
      for(p=1;p<low_num[r];p++)
       {
         low_len[r,p]=low_len[r,p-1] + sqrt( (lastx-low_x[r,p])**2 + (lasty-low_y[r,p])**2 + (lastz-low_z[r,p])**2 );
         lastx=low_x[r,p];  lasty=low_y[r,p];  lastz=low_z[r,p];
       }
      low_totlen[r]=low_len[r,p-1];
    }
   keel_len[0]=0;
   lastx=keel_x[0];  lasty=keel_y[0];  lastz=keel_z[0];
   for(p=1;p<keel_num;p++)
    {
      keel_len[p]=keel_len[p-1] + sqrt( (lastx-keel_x[p])**2 + (lasty-keel_y[p])**2 + (lastz-keel_z[p])**2 );
      lastx=keel_x[p];  lasty=keel_y[p];  lastz=keel_z[p];
    }
   keel_totlen=keel_len[p-1];

   keel_upplen=420*scale_abs; # how far down from the tip do we want the upper hull section to meet the keel?

   fractiondt=1/128;  # step size for upper, covering 0 to <=1 . ie, 0.5 means beginning, middle and end of upper section of rib
   #fractiondt=0.25;  # for quick testing

   curpoint[0]=0;  # keel point
   for(r=1;r<=num_ribs;r++)
      curpoint[r]=0;
   for(fraction=0;fraction<=1;)
    {
      printf "running for fraction %f\n",fraction;
      # here 0 is the keel
      for(r=0;r<=num_ribs;r++)
         printf "%2d ",curpoint[r];
      printf "\n";

      printf "" > "delme3.txt";
      printf "%f %f %f\n",keel_x[curpoint[0] ], keel_y[curpoint[0] ], keel_z[curpoint[0] ] >> "delme3.txt";
      for(r=1;r<=num_ribs;r++)
         printf "%f %f %f\n",upp_x[r,curpoint[r] ], upp_y[r,curpoint[r] ], upp_z[r,curpoint[r] ] >> "delme3.txt";
      r--;  # repeat last line (and later ignore x>max x), closest approx to s[]=0 in spline
      printf "%f %f %f\n",upp_x[r,curpoint[r] ], upp_y[r,curpoint[r] ], upp_z[r,curpoint[r] ] >> "delme3.txt";
      close("delme3.txt");
      system("gawk -f xspline3stdin.awk < delme3.txt > delme4.txt");

      rn=0;
      while( (getline < "delme4.txt")==1 )
       {
         if( $2>upp_x[r,curpoint[r] ] )  # hack in loop avoidance, for the last point, which should be s[]=0
            break;
         row_x[fraction,rn]=$2;  row_y[fraction,rn]=$3;  row_z[fraction,rn]=$4;
         rn++;

         printf "%f %f %f\n",$2,$3,$4 >> "delme2.txt";
       }

      # end of row always missing from xspline3stdin.awk, add it back to the row
      row_x[fraction,rn]=upp_x[r,curpoint[r] ];  row_y[fraction,rn]=upp_y[r,curpoint[r] ];  row_z[fraction,rn]=upp_z[r,curpoint[r] ];
      rn++;

      row_num[fraction]=rn;
      close("delme4.txt");

      #exit;  # stop here for testing initial pass

      fraction+=fractiondt;

      for(r=1;r<=num_ribs;r++)  # find next point on rib for this new fraction
         while( upp_len[r,curpoint[r]] < upp_totlen[r]*fraction && curpoint[r]<upp_num[r] )
            curpoint[r]++;
      lastkeelpoint=curpoint[0];  # save for lower ribs starting point
      while( keel_len[curpoint[0]] < keel_upplen*fraction && curpoint[0]<keel_num )
         curpoint[0]++;
    }

   # and save as obj, storing all the upper rib vertices
   printf "# upper rib vertices\n" >> fn;
   stepdx=1*scale_abs;  # how much to step x by, as the coords in row_x[] are probably higher res than needed
   for(fraction=0;fraction<=1;)
    {
      rn=0;
      if( row_y[fraction,rn] <0 ) row_y[fraction,rn]=0;  # fix hull meeting its mirror overlapping and messing up convex hull property.
      printf "v %f %f %f\n",row_x[fraction,rn], row_y[fraction,rn], row_z[fraction,rn] >> fn;  # original
      row_v[fraction,rn]=vertex++;
      printf "v %f %f %f\n",row_x[fraction,rn], -row_y[fraction,rn], row_z[fraction,rn] >> fn;  # mirror
      vertex++;

      while( rn<row_num[fraction] )
       {
         for(nextrn=rn+1;nextrn<row_num[fraction] && row_x[fraction,rn]+stepdx>row_x[fraction,nextrn];nextrn++)
          ;
         rn=nextrn;
         if( rn<row_num[fraction] )
          {
            if( row_y[fraction,rn] <0 ) row_y[fraction,rn]=0;  # fix hull meeting its mirror overlapping and messing up convex hull property.
            printf "v %f %f %f\n",row_x[fraction,rn], row_y[fraction,rn], row_z[fraction,rn] >> fn;  # original
            row_v[fraction,rn]=vertex++;
            printf "v %f %f %f\n",row_x[fraction,rn], -row_y[fraction,rn], row_z[fraction,rn] >> fn;  # mirror
            vertex++;
          }
         #else
         # {
         #   rnt=row_num[fraction];
         #   printf "v %f %f %f\n",row_x[fraction,rnt], row_y[fraction,rnt], row_z[fraction,rnt] >> fn;
         #   row_v[fraction,rnt]=vertex++;
         # }
       }

      fraction+=fractiondt;
    }

   # lower ribs

   # scale 10 points to fit on the curve of the bow
   keelsteps[int(0)]=0 * scale_abs;
   keelsteps[int(1)]=25 * scale_abs;
   keelsteps[int(2)]=50 * scale_abs;
   keelsteps[int(3)]=85 * scale_abs;
   keelsteps[int(4)]=120 * scale_abs;
   keelsteps[int(5)]=165 * scale_abs;
   keelsteps[int(6)]=210 * scale_abs;
   keelsteps[int(7)]=260 * scale_abs; 
   keelsteps[int(8)]=320 * scale_abs;  
   keelsteps[int(9)]=400 * scale_abs;
   keelsteps[int(10)]=450 * scale_abs;

   # first, copy the last of the upper coards to the first of the lower coards, because some of the upper ribs
   # have been moved slightly, so the upper and lower ribs wont meet.
   for(r=1;r<=num_ribs;r++)
    {
      low_x[r,0]=upp_x[r,upp_num[r]-1];   low_y[r,0]=upp_y[r,upp_num[r]-1];  low_z[r,0]=upp_z[r,upp_num[r]-1];
    }

   # se keel point to last valid step, it has ++ed since it draw a line
   curpoint[0]=lastkeelpoint;

   for(r=1;r<=num_ribs;r++)
      curpoint[r]=0;
   for(fraction=0;fraction<=1;)
    {
      printf "running for fraction %f\n",fraction;
      # here 0 is the keel
      for(r=0;r<=num_ribs;r++)
         printf "%2d ",curpoint[r];
      printf "\n";
      

      printf "" > "delme3.txt";
      printf "%f %f %f\n",keel_x[curpoint[0] ], keel_y[curpoint[0] ], keel_z[curpoint[0] ] >> "delme3.txt";
      for(r=1;r<=num_ribs;r++)
         printf "%f %f %f\n",low_x[r,curpoint[r] ], low_y[r,curpoint[r] ], low_z[r,curpoint[r] ] >> "delme3.txt";
      r--;  # repeat last line (and later ignore x>max x), closest approx to s[]=0 in spline
      printf "%f %f %f\n",low_x[r,curpoint[r] ], low_y[r,curpoint[r] ], low_z[r,curpoint[r] ] >> "delme3.txt";
      close("delme3.txt");
      system("gawk -f xspline3stdin.awk < delme3.txt > delme4.txt");
      rn=0;
      while( (getline < "delme4.txt")==1 )
       {
         if( $2>low_x[r,curpoint[r] ] )  # hack in loop avoidance, for the last point, which should be s[]=0
            break;
         row2_x[fraction,rn]=$2;  row2_y[fraction,rn]=$3;  row2_z[fraction,rn]=$4;
         rn++;
         printf "%f %f %f\n",$2,$3,$4 >> "delme2.txt";
       }

      # end of row always missing from xspline3stdin.awk, add it back to the row
      row2_x[fraction,rn]=low_x[r,curpoint[r] ];  row2_y[fraction,rn]=low_y[r,curpoint[r] ];  row2_z[fraction,rn]=low_z[r,curpoint[r] ];
      rn++;

      row2_num[fraction]=rn;
      close("delme4.txt");

      #exit;  # stop here for testing initial pass

      fraction+=fractiondt;  

      # advance rib points
      for(r=1;r<=num_ribs;r++)
         while( low_len[r,curpoint[r]] < low_totlen[r]*fraction && curpoint[r]<low_num[r] )
             curpoint[r]++;
      # advance keel point
      keelstep=int(fraction*10);
      linearfraction=keelsteps[keelstep] + (keelsteps[keelstep+1]-keelsteps[keelstep])*10*(fraction-keelstep/10);
      printf "keelstep=%f\n",keelstep;
      while( keel_len[curpoint[0]]-keel_upplen < linearfraction && curpoint[0]<keel_num )
         curpoint[0]++;
    }

   # and save as obj, storing all the lower rib vertices
   printf "# lower rib vertices\n" >> fn;
   stepdx=1;  # how much to step x by, as the coords in row_x[] are probably higher res than needed
   for(fraction=0;fraction<=1;)
    {
      rn=0;
      if( row2_y[fraction,rn] <0 ) row2_y[fraction,rn]=0;  # fix hull meeting its mirror overlapping and messing up convex hull property.
      printf "v %f %f %f\n",row2_x[fraction,rn], row2_y[fraction,rn], row2_z[fraction,rn] >> fn;  # original
      row2_v[fraction,rn]=vertex++;
      printf "v %f %f %f\n",row2_x[fraction,rn], -row2_y[fraction,rn], row2_z[fraction,rn] >> fn;  # mirror
      vertex++;

      while( rn<row2_num[fraction] )
       {
         for(nextrn=rn+1;nextrn<row2_num[fraction] && row2_x[fraction,rn]+stepdx>row2_x[fraction,nextrn];nextrn++)
          ;
         rn=nextrn;
         if( rn<row2_num[fraction] )
          {
            if( row2_y[fraction,rn] <0 ) row2_y[fraction,rn]=0;  # fix hull meeting its mirror overlapping and messing up convex hull property.
            printf "v %f %f %f\n",row2_x[fraction,rn], row2_y[fraction,rn], row2_z[fraction,rn] >> fn;  # original
            row2_v[fraction,rn]=vertex++;
            printf "v %f %f %f\n",row2_x[fraction,rn], -row2_y[fraction,rn], row2_z[fraction,rn] >> fn;  # mirror
            vertex++;
          }
       }

      fraction+=fractiondt;
    }

   # vertex in centre of stern, 10mm below the top - so that we can also add a triangle to join the stern faces to the full width
   # faces of the hull top
   printf "v %f %f %f\n", upp_x[num_ribs,0], 0, upp_z[num_ribs,0]-10 >> fn;
   stern_vertex=vertex++;


   printf "# upper rib faces\n" >> fn;
   printf "s off\n" >> fn;  # shading off
   for(fraction=0;fraction<=1-fractiondt;)
    {
      rn=0;
      rn2=0;

      fraction2=fraction+fractiondt;

      while( rn<row_num[fraction] && rn2<row_num[fraction2] )  # move along each fractional row, until the end (the stern)
       {
         if( row_x[fraction,rn] < row_x[fraction2,rn2] )  # move along whichever is more behind
          {
            for(nextrn=rn+1;nextrn<row_num[fraction] && row_x[fraction,rn]+stepdx>row_x[fraction,nextrn];nextrn++)
             ;
            if( nextrn<row_num[fraction] )
             {  # original face and mirror
               printf "f %d %d %d\n", row_v[fraction,rn], row_v[fraction,nextrn], row_v[fraction2,rn2] >> fn;
               printf "f %d %d %d\n", row_v[fraction,rn]+1, row_v[fraction2,rn2]+1, row_v[fraction,nextrn]+1 >> fn;
             }
            rn=nextrn;
          }
         else
          {
            for(nextrn=rn2+1;nextrn<row_num[fraction2] && row_x[fraction2,rn2]+stepdx>row_x[fraction2,nextrn];nextrn++)
             ;
            if( nextrn<row_num[fraction2] )
             {  # original face and mirror
               printf "f %d %d %d\n", row_v[fraction2,nextrn], row_v[fraction2,rn2], row_v[fraction,rn] >> fn;
               printf "f %d %d %d\n", row_v[fraction2,nextrn]+1, row_v[fraction,rn]+1, row_v[fraction2,rn2]+1 >> fn;
             }
            rn2=nextrn;
          }

       }

      fraction+=fractiondt;
      fraction2+=fractiondt;
    }

   # here the 2 refers to the lower part of the rib, except for rn2 - which is the same use as above
   printf "# lower rib faces\n" >> fn;
   for(fraction=0;fraction<=1-fractiondt;)
    {
      rn=0;
      rn2=0;

      fraction2=fraction+fractiondt;

      while( rn<row2_num[fraction] && rn2<row2_num[fraction2] )  # move along each fractional row, until the end (the stern)
       {
         if( row2_x[fraction,rn] < row2_x[fraction2,rn2] )  # move along whichever is more behind
          {
            for(nextrn=rn+1;nextrn<row2_num[fraction] && row2_x[fraction,rn]+stepdx>row2_x[fraction,nextrn];nextrn++)
             ;
            if( nextrn<row2_num[fraction] )
             {  # original face and mirror
               printf "f %d %d %d\n", row2_v[fraction,rn], row2_v[fraction,nextrn], row2_v[fraction2,rn2] >> fn;
               printf "f %d %d %d\n", row2_v[fraction,rn]+1, row2_v[fraction2,rn2]+1, row2_v[fraction,nextrn]+1 >> fn;
             }
            rn=nextrn;
          }
         else
          {
            for(nextrn=rn2+1;nextrn<row2_num[fraction2] && row2_x[fraction2,rn2]+stepdx>row2_x[fraction2,nextrn];nextrn++)
             ;
            if( nextrn<row2_num[fraction2] )
             {  # original face and mirror
               printf "f %d %d %d\n", row2_v[fraction2,nextrn], row2_v[fraction2,rn2], row2_v[fraction,rn] >> fn;
               printf "f %d %d %d\n", row2_v[fraction2,nextrn]+1, row2_v[fraction,rn]+1, row2_v[fraction2,rn2]+1 >> fn;
             }
            rn2=nextrn;
          }

       }

      fraction+=fractiondt;
      fraction2+=fractiondt;
    }


   # save a stern, stern_vertex joined with all pairs of upp_v and then low_v
   printf "# stern faces\n" >> fn;

   printf "f %d %d %d\n", stern_vertex, row_v[0,row_num[0]-1], row_v[0,row_num[0]-1]+1 >> fn;  # face above stern_vertex

   for(fraction=0;fraction<=1-fractiondt;fraction+=fractiondt)  # upper ribs
    {
      fraction2=fraction+fractiondt;
      rn=row_num[fraction]-1;
      rn2=row_num[fraction2]-1;
      
      printf "f %d %d %d\n", stern_vertex, row_v[fraction2,rn2], row_v[fraction,rn] >> fn;
      printf "f %d %d %d\n", stern_vertex, row_v[fraction,rn]+1, row_v[fraction2,rn2]+1 >> fn;
      lastv=row_v[fraction2,rn2];
    }

   printf "# stern faces, joining the mirrored half\n" >> fn;
   printf "f %d %d %d\n", stern_vertex, lastv, row2_v[0,row2_num[0]-1] >> fn;
   printf "f %d %d %d\n", stern_vertex, lastv+1, row2_v[0,row2_num[0]-1]+1 >> fn;
 
   for(fraction=0;fraction<=1-fractiondt;fraction+=fractiondt)  # lower ribs
    {
      fraction2=fraction+fractiondt;
      rn=row2_num[fraction]-1;
      rn2=row2_num[fraction2]-1;
      
      printf "f %d %d %d\n", stern_vertex, row2_v[fraction2,rn2], row2_v[fraction,rn] >> fn;
      printf "f %d %d %d\n", stern_vertex, row2_v[fraction,rn]+1, row2_v[fraction2,rn2]+1 >> fn;  # mirror
    }

   # doesnt meet in the middle, so add another face that joins the original and the mirror 
   printf "f %d %d %d\n", stern_vertex, row2_v[fraction2,rn2], row2_v[fraction2,rn2]+1 >> fn;


   # and now add a lid to close the object
   printf "# lid\n" >> fn;
   fraction=0;  # top row
   rn=0;
   rn2=0;

   while( rn<row_num[fraction] && rn2<row_num[fraction] )
    {
      if( row_x[fraction,rn] < row_x[fraction,rn2] )  # move along whichever is more behind
       {
         for(nextrn=rn+1;nextrn<row_num[fraction] && row_x[fraction,rn]+stepdx>row_x[fraction,nextrn];nextrn++)
          ;
         if( nextrn<row_num[fraction] )
            printf "f %d %d %d\n", row_v[fraction,nextrn], row_v[fraction,rn], row_v[fraction,rn2]+1 >> fn;
         rn=nextrn;
       }
      else
       {
         for(nextrn=rn2+1;nextrn<row_num[fraction] && row_x[fraction,rn2]+stepdx>row_x[fraction,nextrn];nextrn++)
          ;
         if( nextrn<row_num[fraction] )
            printf "f %d %d %d\n", row_v[fraction,rn]+1, row_v[fraction,nextrn]+1, row_v[fraction,rn] >> fn;
         rn2=nextrn;
       }
    }

   # and then add lower to mirrored lower rib faces, as meshlab just cant make the object
   # water tight without creating other artifacts (like changing the shape...).
   printf "# left to right keel joiner\n" >> fn;
   fraction=1;  # bottow row
   rn=0;
   rn2=0;

   while( rn<row2_num[fraction] && rn2<row2_num[fraction] )
    {
      if( row2_x[fraction,rn] < row2_x[fraction,rn2] )  # move along whichever is more behind
       {
         for(nextrn=rn+1;nextrn<row2_num[fraction] && row2_x[fraction,rn]+stepdx>row2_x[fraction,nextrn];nextrn++)
          ;
         if( nextrn<row2_num[fraction] )
          {
            #if( row2_y[fraction,nextrn] > row2_y[fraction,rn2] )  # this is never true
            #   printf "f %d %d %d\n", row2_v[fraction,nextrn], row2_v[fraction,rn], row2_v[fraction,rn2]+1 >> fn;
            #else
               printf "f %d %d %d\n", row2_v[fraction,nextrn], row2_v[fraction,rn2]+1,row2_v[fraction,rn] >> fn;
          }
         rn=nextrn;
       }
      else
       {
         for(nextrn=rn2+1;nextrn<row2_num[fraction] && row2_x[fraction,rn2]+stepdx>row2_x[fraction,nextrn];nextrn++)
          ;
         if( nextrn<row2_num[fraction] ) 
          {
            #if( row2_y[fraction,nextrn] > row2_y[fraction,rn] )
               printf "f %d %d %d\n", row2_v[fraction,rn2]+1, row2_v[fraction,rn], row2_v[fraction,nextrn]+1 >> fn;
            #else
            #   printf "f %d %d %d\n", row2_v[fraction,rn2]+1, row2_v[fraction,nextrn]+1, row2_v[fraction,rn] >> fn;
          }
         rn2=nextrn;
       }
    }



   close(fn);


   exit;

   minz=maxz=keel_z[1];
   for(p in keel_z)
      if( keel_z[p]<minz )minz=keel_z[p];
      else if( keel_z[p]>max_z )maxz=keel_z[p];
   for(r=1;r<=num_ribs;r++)
      for(p=0;p<upp_num[r];p++)
         if( upp_z[r,p]<minz )minz=upp_z[r,p];
         else if( upp_z[r,p]>max_z )maxz=upp_z[r,p];

   printf "minz=%f, maxz=%f\n",minz,maxz;



   for(z=maxz;z>=minz;z-=1)
    {
      sol=0;
      for(p=1;p<keel_num;p++)
         if( keel_z[p-1]<=z && keel_z[p]>=z )
          {
            y=inter( keel_y[p-1], keel_z[p-1],  keel_y[p], keel_z[p], z );  # always 0 for the keel
            x=inter( keel_x[p-1], keel_z[p-1],  keel_x[p], keel_z[p], z );
            sol=1;
            break;
          }
      if( sol==0 )
         printf "no solution for z=%f\n",z;

      printf "z=%f, p=%d, (%0.3f,%0.3f,%0.3f)-(%0.3f,%0.3f,%0.3f) => (%0.3f,%0.3f)\n",
         z,p, keel_x[p-1], keel_y[p-1], keel_z[p-1],  keel_x[p], keel_y[p], keel_z[p],   x,y;

      for(rib=1;rib<=num_ribs && sol>0;rib++)
       {
       }
    }

 }' delme.txt