!**********************************************

  BASIC GEOMETRY MODULE caldspmnt(
   REF hull_id >"@t16 Select hull surface !";
   REF wp_id   >"@t16 Select waterplane !";
   INT ns:=10  >"Enter number of double Simpson slices !");

!*   Numerical calculation of displacement,
!*   center of gravity, waterplane area and
!*   moments of inertia.
!*
!*   Calculates and stores:
!*
!*   LWL(2)          = Length in waterline
!*   DISPLACEMENT(4) = Volume
!*   WPAREA(6)       = Waterplane area
!*
!*   See file "Datatypes" for further information.
!*
!*   (C)Microform AB, J.Kjellander
!*      Stephan Helma
!*
!**********************************************

 INT    i,status;
 FLOAT  dx,lwl,ywl,
        x1,x2,x3,
        area1,area2,area3,
	beam1,beam2,beam3,
	slice_volume,total_volume,fx,
	slice_wparea,total_wparea,bx,bxx,bxxx,bbb,
        k(6);
 VECTOR p1,p2,tmp;
 REF    stern_id;
 STRING s*132;

BEGINMODULE

!*
!***Check number of Simpson slices.
!*
   if ns < 1 then
     exit("Number of double Simpson slices must be >0 !");
   endif;
!*
!***Global ID of hull surface.
!*
   hull_id:=global_ref(hull_id,9);
!*
!***Get y-value of waterplane.
!*
   ywl:=getflt(wp_id,2);
!*
!***ID for waterplane b_plane.
!*
   s:=rstr(wp_id);
   wp_id:=rval(s+"#1");  
!*
!***Stern surface ID.
!*
   s:=rstr(hull_id);
   stern_id:=rval(s+"#50");  
!*
!***Create actual waterplane curve by intersecting
!***hullsurface with waterplane surface b_plane.
!*
   cur_int(#1,hull_id,wp_id:BLANK=0,PEN=5);
!*
!***Where does waterline begin and end.
!*
   p1:=startp(#1);
   p2:=endp(#1);
   lin_free(#2,p1,p2:PEN=5);

   if p2.x < p1.x then
     tmp:=p1;
     p1:=p2;
     p2:=tmp;
   endif;

   lwl:=p2.x-p1.x;
!*
!***Calculate dx, so that a multiple of
!***it fits into LWL. If dx < 1 correct ns
!***and store it into database.
!*
   dx:=lwl/(2.0*ns);
   if dx < 1 then
     ns:=trunc(lwl/2.0);
     dx:=lwl/(2.0*ns);
   endif;
!*
!***Slice the hull along the waterline.
!***fx  is the volume of each slice
!***    multiplied with x-coordinate.
!***bx  the beam multiplied with x-coord.
!***bxx the beam multiplied with x-coord^2.
!***bbb the beam^3
!*
   total_volume:=0.0;
   fx:=0.0;
   total_wparea:=0.0;
   bx:=0.0;
   bxx:=0.0;
   bbb:=0.0;

   area3:=0.0;
   beam3:=0.0;

   for i:=1 to ns do
!     x:=p1.x+i*2.0*dx;
     x1:=p1.x+(i-1)*2.0*dx;
     x2:=x1+dx;
     x3:=x1+2.0*dx;
     area1:=area3;
     beam1:=beam3;
     
     part(#3,sectarea(hull_id,wp_id,ywl,x2,beam2,area2));
     if i = ns then
       area3:=0.0;
       beam3:=0.0;
     else
       part(#4,sectarea(hull_id,wp_id,ywl,x3,beam3,area3));
     endif;
!FIXME: Is it possible to find out if sectarea cannot finish
!       successfully, because there is no intersection? That
!       would even work for a plan stern (or bow, like the
!       Optimist.

!*
!***Calculate enclosed area. We use Simpson's rule to get
!***a better result. We are at position x with area1. area2
!***is the area of the middle section and area3 is the area
!***of the last section for Simpson's rule.
!***The first area of the first slice must be zero as the
!***third area of the last slice. (IF THE HULL DOESN'T END ABRUPTLY)
!*

     slice_volume:=area1+4.0*area2+area3;
     total_volume:=total_volume+slice_volume;
!     fx:=fx+slice_volume*(x-dx);
!     fx:=fx+(area1*(x-2.0*dx)+4.0*area2*(x-dx)+area3*x);
     fx:=fx+(area1*x1+4.0*area2*x2+area3*x3);
     slice_wparea:=beam1+4.0*beam2+beam3;
     total_wparea:=total_wparea+slice_wparea;
     bx:=bx+(beam1*x1+4.0*beam2*x2+beam3*x3);
!     bx:=bx+(beam1*(x-2.0*dx)+4.0*beam2*(x-dx)+beam3*x);
     bxx:=bxx+(beam1*x1*x1+4.0*beam2*x2*x2+beam3*x3*x3);
!     bxx:=bxx+(beam1*(x-2.0*dx)**2+4.0*beam2*(x-dx)**2+beam3*x**2);
!FIXME: for the moment we assume symmetrical waterplane
!       (=upright condition)!!
!     bbb:=beam1**3+4.0*beam2**3+beam3**3;
     bbb:=bbb+((beam1/2.0)**3+4.0*(beam2/2.0)**3+(beam3/2.0)**3);

     pop_pmt();
     psh_pmt("Volume of slice no:"+str(i,-1,0)+" is "+
              str(slice_volume/1E9,-1,3)+" m³");

   endfor;

   pop_pmt();
   psh_pmt("Storing results in Database !");
   total_volume:=total_volume*dx/3.0;
   fx:=fx*dx/3.0;
   total_wparea:=total_wparea*dx/3.0;
   bx:=bx*dx/3.0;
   bxx:=bxx*dx/3.0;
   bbb:=bbb*dx/3.0/3.0;   
   
   k(1):=total_volume;
   k(2):=fx/total_volume;
   k(3):=99999;
   k(4):=99999;
   status:=deldat_gm("DISPLACEMENT");
   status:=putdat_gm("DISPLACEMENT",4,k);
   k(1):=total_wparea;
   k(2):=bx/total_wparea;
   k(3):=ywl;
   k(4):=0.0;
   k(5):=bxx-bx*bx/total_wparea;
   k(6):=bbb;
   status:=deldat_gm("WPAREA");
   status:=putdat_gm("WPAREA",6,k);

   pop_pmt();

ENDMODULE

!*****************************************************