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

  BASIC GEOMETRY MODULE hydrstat(
   REF   hull_id      >"@t16 Select hull surface !";
   FLOAT design_displ >"Enter design displacement in m³ !";
   FLOAT heel:=0.0    >"Heel angle !";
!!   FLOAT dx:=1000   >"Step size !");
   INT   ns:=10       >"Enter number of double Simpson slices !");

!*   Balance the hull through numerical calculation
!*   of displacement and calculate center of
!*   gravity, waterplane area and moments.
!*
!*   Calculates and stores:
!*
!*   not up-to-date! See file "Datatypes".
!*   LWL(2)            = Length of waterline
!*   SIMPSON_SLICES(1) = Number of double Simpson slices
!*
!*   See file "Datatypes" for further information.
!*
!*   (C)Microform AB, J.Kjellander
!*      Stephan Helma
!*
!**********************************************

!! INT    i,ns,font,ns_geo,ns_gra,status;
 INT    i,status;
 FLOAT  ywl,ywl_new,
        wparea_new,
	displ_new,
	xtolerance,
        k(6);
 REF    wp_id,stern_id;
 STRING s*132;


BEGINMODULE
pop_pmt();
!*
!***Check number of Simpson slices.
!*
   if ns < 1 then
     exit("Number of double Simpson slices must be >0 !");
   endif;
   k(1):=ns;
   status:=deldat_gm("SIMPSON_SLICES");
   status:=putdat_gm("SIMPSON_SLICES",1,k);

!*
!***Global ID of hull surface.
!*
   hull_id:=global_ref(hull_id,9);
!*
!***Stern surface ID.
!*
   s:=rstr(hull_id);
   stern_id:=rval(s+"#50");  

!*
!***Get vertical position of actual waterplane
!***as a starting point for balancing.
!*
   status:=getdat_gm("LWL",k);
   if status < 2 then
     exit("No waterplane found in database !");
   endif;
   ywl:=k(2);
   design_displ:=design_displ*1E9;
   xtolerance:=1.0;
   i:=1;

!*
!***Now balance the hull. We use Newton's
!***method to speed up balancing. Since the
!***area of waterplane is the first
!***derivative of the displacement it's quite
!***easy. If the change in ywl is less than
!***XTOLERANCE we got the floating WL.
!*

start_balancing:
   psh_pmt("Balancing the hull: step "+str(i,-1,0)+" !");
   part(#1,waterpl(hull_id,ywl+100,heel):BLANK=1);
   wp_id:=global_ref(#1);
   part(#2,caldspmnt(hull_id,wp_id,ns));
   status:=getdat_gm("DISPLACEMENT",k);
   displ_new:=k(1);
   status:=getdat_gm("WPAREA",k);
   wparea_new:=k(1);
   ywl_new:=ywl+(design_displ-displ_new)/wparea_new;
   pop_pmt();
   if abs(ywl-ywl_new) > xtolerance then
     ywl:=ywl_new;
     del(#1);
     del(#2);
     i:=i+1;
     goto start_balancing;
   endif;

!*
!***Max beam and area of section.
!*
   part(#5,calmsect(hull_id,ywl));

!*
!***Hull + stern area.
!*
   part(#6,calwsurf(hull_id));

ENDMODULE

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