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

  BASIC GEOMETRY MODULE dispmnt(
   REF hull_id >"@t16 Select hull surface !";
   FLOAT dx:=1000  >"Step size !");

!*   Numerical calculation of displacement and
!*   center of gravity.
!*
!*   Calculates and stores:
!*
!*   LWL(1)          = Length in waterline in meters
!*   DISPLACEMENT(1) = Volume im m3
!*   DISPLACEMENT(2) = CB X-station in meters
!*   CP(1)           = Prismatic coefficient
!*   DLRATIO(1)      = Displacement/LWL ratio in pounds/feet
!*   TBCOMFORT(1)    = Ted Brewers Comfort factor
!*   DISPSPEED(1)    = Max theoretical displacement speed in knots
!*   HULLAREA(1)     = Hull area exkl. stern in m2
!*   HULLAREA(2)     = Stern area in m2
!*
!*   (C)Microform AB, J.Kjellander
!*
!**********************************************

 INT    i,ns,font,ns_geo,ns_gra,status;
 FLOAT  x,ywl,area1,area2,slice_area,lwl,speed,
        mb_area,slice_volume,total_volume,fx,dl,comfort,
        k(2);
 VECTOR p1,p2,pwl,pa,pb,tmp;
 REF    wp_id,stern_id;
 STRING s*132;

string sjk*1; int njk;

BEGINMODULE

!*
!***Check stepsize.
!*
   if dx < 1 then
     exit("Step size is too small !");
   endif;
!*
!***Global ID of hull surface.
!*
   hull_id:=global_ref(hull_id,9);
!*
!***Global ID of waterplane.
!*
   set_root_gm();
   wp_id:=get_next_gm(1,16384,"waterpl");
   if wp_id = #0.1 then
     exit("No waterplane found !");
   endif;
!*
!***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;
!*
!***Some text output.
!*
   ns:=trunc(lwl/dx)+1;
!*
!***Slice the hull along the waterline. fx is the volume of each
!***slice multiplied with X-coordinate.
!*
   fx:=0.0;
   total_volume:=0.0;
   i:=1;
   x:=p1.x+dx;
loop: 
   part(#3,sectarea(hull_id,wp_id,ywl,x,area2));
!*
!***Calculate enclosed area. area1 is the area of the
!***previous section and area2 is the area of this section.
!***Full area is twice as big as we only operate on
!***the left half of the hull. The area of the current
!***slice is approximated by taking the mean value of
!***area1 and area2. Special treatment is done for the
!***first slice to increase accuracy. It's area is
!***calculated by assuming a previous section with
!***area = 0. Better approximations could be invented
!***but this will do for now.
!*
   if i = 1 then
     area1:=0;
   else
     area1:=area2;
   endif;

   slice_area:=(area1+area2)/2.0;
   slice_volume:=slice_area*dx;
   total_volume:=total_volume+slice_volume;
   fx:=fx+slice_volume*(x-dx/2);

   pop_pmt();
   psh_pmt("Volume of slice no:"+str(i,-1,0)+" is "+
            str(slice_volume/1E6,-1,0)+" kg");
!*
!***Next section.
!*
   if x+dx < p2.x then
     x:=x+dx;
     i:=i+1;
     goto loop;
   endif;
!*
!***Last slice is not in the loop above. We use the
!***same approximation for the last slice as for the
!***first slice exept that the thickness of the last
!***slice has to be calculated independantly since it
!***is usually < dx.
!*
   dx:=p2.x-x;
   slice_area:=area2/2.0;
   slice_volume:=slice_area*dx;
   total_volume:=total_volume+slice_volume;
   fx:=fx+slice_volume*(x+dx/2);

   pop_pmt();
   psh_pmt("Volume of slice no:"+str(i+1,-1,0)+" is "+
            str(slice_volume/1E6,-1,0)+" kg");

   k(1):=total_volume/1E6;
   k(2):=fx/total_volume/1E3;
   status:=deldat_gm("DISPLACEMENT");
   status:=putdat_gm("DISPLACEMENT",2,k);
!*
!***To calculate prismatic coefficient we need the max
!***beam of the waterline so we can calculate the area
!***of the max beam section. NOTE, this is not neccesary
!***the same station as MAXBEAM of the entire hull !
!*
   getcurh(#1,font,dl,ns_geo,ns_gra);
   pb:=vec(0,0,0);

   for i:=0 to 25*ns_geo do
     pa:=on(#1,i/25);
     if pa.z > pb.z then
       pb:=pa;
     endif;
   endfor;

   part(#4,sectarea(hull_id,wp_id,ywl,pb.x,mb_area));
   pop_pmt();

   k(1):=total_volume/(mb_area*lwl);

   status:=deldat_gm("CP");
   status:=putdat_gm("CP",1,k);
!*
!***D/L ratio.
!*
   dl:=total_volume/1E6/0.4536/2240/((0.01*lwl/1E3/0.3048)**3);

   k(1):=dl;
   status:=deldat_gm("DLRATIO");
   status:=putdat_gm("DLRATIO",1,k);
!*
!***Max theoretical displacement speed.
!*
   speed:=1.34*sqrt(lwl/1E3/0.3048);

   k(1):=speed;
   status:=deldat_gm("DISPSPEED");
   status:=putdat_gm("DISPSPEED",1,k);
!*
!***Ted Brewers "comfort factor". To that we nead LOA as well.
!***Note that there is a typing error on page 8 in his book.
!***65 should be 0.65.
!*
   comfort:=total_volume/1E6/0.4536/
            (0.65*(0.7*lwl/1E3/0.3048+0.3*53)*(2*pb.z/1E3/0.3048)**1.333);

   k(1):=comfort;
   status:=deldat_gm("TBCOMFORT");
   status:=putdat_gm("TBCOMFORT",1,k);
!*
!***Hull + stern area.
!*
   psh_pmt("Calculating hull + stern area !");
   k(1):=(2*surface_area(hull_id))/1E6;
   k(2):=0; ! surface_area(stern_id)/1E6;
   status:=deldat_gm("HULLAREA");
   status:=putdat_gm("HULLAREA",2,k);
   pop_pmt();

ENDMODULE

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