From 718d422b60cd29b53c9404c641c88d3f553843b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Alejandro=20Iv=C3=A1=C3=B1ez=20Encinas?= Date: Tue, 17 Oct 2017 19:59:44 +0200 Subject: [PATCH 01/20] Changed Stability & added block case I've changed the two variables that were missplaced and I've added blockCase_compare.html. I've calculated a lot of parameters of the blockCase and printed them in a table. I'll add the rest. --- build/Vessel.js | 6 +- examples/blockCase_compare.html | 371 ++++++++++++++++++ .../data/ship_specifications/blockCase.json | 67 ++++ source/CoreClasses/Ship.js | 6 +- 4 files changed, 444 insertions(+), 6 deletions(-) create mode 100644 examples/blockCase_compare.html create mode 100644 examples/data/ship_specifications/blockCase.json diff --git a/build/Vessel.js b/build/Vessel.js index 0501576..bbc1386 100644 --- a/build/Vessel.js +++ b/build/Vessel.js @@ -701,9 +701,9 @@ Object.assign(Ship.prototype, { vol = Lwl * B * T * cb; } let KG = this.getWeight(shipState).cg.z; - let I = ha.Iywp * 1000; - let KB = 0.52 * T; - let BM = I / vol; + let I = ha.Iywp; + let BM = 0.52 * T; + let KB = I / vol; let GM = KB + BM - KG; return {GM, KB, BM, KG}; } diff --git a/examples/blockCase_compare.html b/examples/blockCase_compare.html new file mode 100644 index 0000000..00315c3 --- /dev/null +++ b/examples/blockCase_compare.html @@ -0,0 +1,371 @@ + + + + + + + + + + + + + Ship visualization with specification + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + +
+ + + + +
+

Ship as object

+
+ +
+ + + +

This app asks for a JSON file + and gives back the pretty form of the JSON element. It also plots a 3D image of the ship. Sample data is loaded at startup.

+ +

Developed by: Elias Hasle, Vicente Alejandro Iváñez Encinas and Henrique M. Gaspar. Visualization made using three.js.

+ +
+ +
+

Input

+ +

Contents of the file:

+ 

+                

+            

+            
+            

+                
3D orbit view of hull and parts

+			

+			

+            

+            
+            

+                
Calculations

+                
GM : m.

+                
+			

+            

+                
Comparation

+                
+                    
+                        
+                            
+                            
+                                                  
+                                             
+                        
+                    
+                    
+                    
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                        
+                        
+                    
+                
+                
VariableHand valueCalculated value% diff
Displacement (t)
Waterline Lenght (m)
Maximum Waterline Breath (m)
Wetted area (m2)
Waterplane Area (m2)
Cb
LCBx (m)
LCFx (m)
KB (m)
KG (m)
BM (m)
GM (m)
+                           

+			

+        

+    

+        
+    
+    
+    
+    
+    
+    
+      
+
+
+
+
\ No newline at end of file
diff --git a/examples/data/ship_specifications/blockCase.json b/examples/data/ship_specifications/blockCase.json
new file mode 100644
index 0000000..b3a4601
--- /dev/null
+++ b/examples/data/ship_specifications/blockCase.json
@@ -0,0 +1,67 @@
+{
+"attributes": {
+},
+"designState": {
+	"calculationParameters": {
+		"LWL_design": 20,
+		"Draft_design": 2,
+		"Cb_design": 1,
+		"speed": 12,
+	    "crew" : 20,
+		"K" : 0.032,
+		"MCR" : 10000,
+		"SFR" : 0.000215,
+		"Co" : 0.3, 
+		"tripDuration": 40
+	},
+	"objectOverrides": {
+		"common": {
+			"fullness": 0.7
+		}
+	}
+},
+"structure": {
+	"hull": {
+		"attributes": {
+			"LOA": 20,
+			"BOA": 10,
+			"Depth": 4,
+			"APP": 0
+		},
+		"halfBreadths": {
+			"waterlines": [0, 1],
+			"stations": [0, 1],
+			"table": [[0, 0], [1, 1]]
+		},
+		"buttockHeights": {}
+	},
+	"decks": {
+		"WheelHouseTop": {
+			"zFloor": 4,
+			"thickness": 0.3,
+			"xAft": 0,
+			"xFwd": 20,
+			"yCentre": 0,
+			"breadth": 10,
+			"density": 1500
+		},
+		"BridgeDeck": {
+			"zFloor": 2,
+			"thickness": 0.3,
+			"xAft": 0,
+			"xFwd": 20,
+			"yCentre": 0,
+			"breadth": 5,
+			"density": 1500
+		}
+		
+	},
+	"bulkheads": {
+	}
+},
+"baseObjects": 
+[{"id":"TankL2.1B6.13H4.3fdundefinedFtank1.stl","affiliations":{},"boxDimensions":{"length":4,"breadth":6,"height":5},"capabilities":{},"file3D":"tank1.stl","baseState":{"fullness":1},"weightInformation":{"contentDensity":1000,"volumeCapacity":120,"lightweight":553.539,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}},{"id":"TankL2.1B6.13H4.3fdundefinedFundefined","affiliations":{},"boxDimensions":{"length":5,"breadth":4,"height":2},"capabilities":{},"baseState":{"fullness":0.5},"weightInformation":{"contentDensity":1000,"volumeCapacity":20,"lightweight":553.539,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}},{"id":"TankL4.199999999999999B6.7H4.3fdundefinedFtank2.stl","affiliations":{},"boxDimensions":{"length":10,"breadth":4,"height":2},"capabilities":{},"file3D":"tank2.stl","baseState":{"fullness":0.5},"weightInformation":{"contentDensity":1000,"volumeCapacity":48,"lightweight":1210.0199999999998,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}}],
+
+"derivedObjects": [{"id":"Tank1","baseObject":"TankL2.1B6.13H4.3fdundefinedFtank1.stl","affiliations":{"Deck":"WheelHouseTop","SFI":"101"},"referenceState":{"xCentre":2,"yCentre":0,"zBase":4}},{"id":"Tank2","baseObject":"TankL2.1B6.13H4.3fdundefinedFundefined","affiliations":{"Deck":"BridgeDeck","SFI":"102"},"referenceState":{"xCentre":5,"yCentre":0,"zBase":2}},{"id":"Tank3","baseObject":"TankL4.199999999999999B6.7H4.3fdundefinedFtank2.stl","affiliations":{"Deck":"BridgeDeck","SFI":"103"},"referenceState":{"xCentre":14,"yCentre":0,"zBase":2}}]
+
+}
\ No newline at end of file
diff --git a/source/CoreClasses/Ship.js b/source/CoreClasses/Ship.js
index 1885542..c1d5796 100644
--- a/source/CoreClasses/Ship.js
+++ b/source/CoreClasses/Ship.js
@@ -94,9 +94,9 @@ Object.assign(Ship.prototype, {
             vol = Lwl * B * T * cb;
         }
         let KG = this.getWeight(shipState).cg.z;
-        let I = ha.Iywp * 1000;
-        let KB = 0.52 * T;
-        let BM = I / vol;
+        let I = ha.Iywp;
+        let BM = 0.52 * T;
+        let KB = I / vol;
         let GM = KB + BM - KG;
         return {GM, KB, BM, KG};
 	}

From b5c44aa031628657851f7873253016a61da58898 Mon Sep 17 00:00:00 2001
From: EliasHasle 
Date: Wed, 18 Oct 2017 09:10:31 +0200
Subject: [PATCH 02/20] added @MrEranwe 's new files, and some modifications

---
 build/Vessel.js                               |   5 +-
 examples/blockCase_compare.html               | 371 ++++++++++++++++++
 .../data/ship_specifications/blockCase.json   |  67 ++++
 source/CoreClasses/Ship.js                    |   3 +-
 4 files changed, 443 insertions(+), 3 deletions(-)
 create mode 100644 examples/blockCase_compare.html
 create mode 100644 examples/data/ship_specifications/blockCase.json

diff --git a/build/Vessel.js b/build/Vessel.js
index 0501576..29040a9 100644
--- a/build/Vessel.js
+++ b/build/Vessel.js
@@ -1,4 +1,4 @@
-//Vessel.js library, built 2017-10-17 12:23:09.339097, Checksum: 5b20e5739dc8e9c7a5c69904053ad53c
+//Vessel.js library, built 2017-10-18 09:07:47.781961, Checksum: 98f5cf6e023f4efb329738cead04ff77
 /*
 Import like this in HTML:
 
@@ -690,6 +690,7 @@ Object.assign(Ship.prototype, {
 		console.info("Calculated draft: %.2f", t);
 		return t;
 	},
+	//Should separate between longitudinal and transverse GM too
     calculateStability(shipState){
         let T = this.calculateDraft(shipState);
         let ha = this.structure.hull.calculateAttributesAtDraft(T);
@@ -701,7 +702,7 @@ Object.assign(Ship.prototype, {
             vol = Lwl * B * T * cb;
         }
         let KG = this.getWeight(shipState).cg.z;
-        let I = ha.Iywp * 1000;
+        let I = ha.Iywp;
         let KB = 0.52 * T;
         let BM = I / vol;
         let GM = KB + BM - KG;
diff --git a/examples/blockCase_compare.html b/examples/blockCase_compare.html
new file mode 100644
index 0000000..00315c3
--- /dev/null
+++ b/examples/blockCase_compare.html
@@ -0,0 +1,371 @@
+
+
+
+
+
+
+    
+    
+    
+    
+    
+
+    Ship visualization with specification
+
+    
+    
+    
+    
+    
+
+    
+     
+	
+	
+	
+	
+	
+    
+	
+    
+
+	
+	
+
+	
+	
+	
+	
+	
+
+    
+    
+    
+
+
+
+
+
+    
+    
+            
+    
+    

+        
+        
+        

+            
+            
+            
+            
+             

+                
Ship as object

+            
  
+            
+			

+                
+                
+                
+                
This app asks for a JSON file 
+                and gives back the pretty form of the JSON element. It also plots a 3D image of the ship. Sample data is loaded at startup.

+				                                    
+				
Developed by: Elias Hasle, Vicente Alejandro Iváñez Encinas and Henrique M. Gaspar. Visualization made using three.js.

+
+            

+			
+            

+                
Input

+                
+                
Contents of the file:

+                

+                

+            

+            
+            

+                
3D orbit view of hull and parts

+			

+			

+            

+            
+            

+                
Calculations

+                
GM : m.

+                
+			

+            

+                
Comparation

+                
+                    
+                        
+                            
+                            
+                                                  
+                                             
+                        
+                    
+                    
+                    
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                            
+                        
+                        
+                            
+                            
+                            
+                        
+                        
+                    
+                
+                
VariableHand valueCalculated value% diff
Displacement (t)
Waterline Lenght (m)
Maximum Waterline Breath (m)
Wetted area (m2)
Waterplane Area (m2)
Cb
LCBx (m)
LCFx (m)
KB (m)
KG (m)
BM (m)
GM (m)
+                           

+			

+        

+    

+        
+    
+    
+    
+    
+    
+    
+      
+
+
+
+
\ No newline at end of file
diff --git a/examples/data/ship_specifications/blockCase.json b/examples/data/ship_specifications/blockCase.json
new file mode 100644
index 0000000..b3a4601
--- /dev/null
+++ b/examples/data/ship_specifications/blockCase.json
@@ -0,0 +1,67 @@
+{
+"attributes": {
+},
+"designState": {
+	"calculationParameters": {
+		"LWL_design": 20,
+		"Draft_design": 2,
+		"Cb_design": 1,
+		"speed": 12,
+	    "crew" : 20,
+		"K" : 0.032,
+		"MCR" : 10000,
+		"SFR" : 0.000215,
+		"Co" : 0.3, 
+		"tripDuration": 40
+	},
+	"objectOverrides": {
+		"common": {
+			"fullness": 0.7
+		}
+	}
+},
+"structure": {
+	"hull": {
+		"attributes": {
+			"LOA": 20,
+			"BOA": 10,
+			"Depth": 4,
+			"APP": 0
+		},
+		"halfBreadths": {
+			"waterlines": [0, 1],
+			"stations": [0, 1],
+			"table": [[0, 0], [1, 1]]
+		},
+		"buttockHeights": {}
+	},
+	"decks": {
+		"WheelHouseTop": {
+			"zFloor": 4,
+			"thickness": 0.3,
+			"xAft": 0,
+			"xFwd": 20,
+			"yCentre": 0,
+			"breadth": 10,
+			"density": 1500
+		},
+		"BridgeDeck": {
+			"zFloor": 2,
+			"thickness": 0.3,
+			"xAft": 0,
+			"xFwd": 20,
+			"yCentre": 0,
+			"breadth": 5,
+			"density": 1500
+		}
+		
+	},
+	"bulkheads": {
+	}
+},
+"baseObjects": 
+[{"id":"TankL2.1B6.13H4.3fdundefinedFtank1.stl","affiliations":{},"boxDimensions":{"length":4,"breadth":6,"height":5},"capabilities":{},"file3D":"tank1.stl","baseState":{"fullness":1},"weightInformation":{"contentDensity":1000,"volumeCapacity":120,"lightweight":553.539,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}},{"id":"TankL2.1B6.13H4.3fdundefinedFundefined","affiliations":{},"boxDimensions":{"length":5,"breadth":4,"height":2},"capabilities":{},"baseState":{"fullness":0.5},"weightInformation":{"contentDensity":1000,"volumeCapacity":20,"lightweight":553.539,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}},{"id":"TankL4.199999999999999B6.7H4.3fdundefinedFtank2.stl","affiliations":{},"boxDimensions":{"length":10,"breadth":4,"height":2},"capabilities":{},"file3D":"tank2.stl","baseState":{"fullness":0.5},"weightInformation":{"contentDensity":1000,"volumeCapacity":48,"lightweight":1210.0199999999998,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}}],
+
+"derivedObjects": [{"id":"Tank1","baseObject":"TankL2.1B6.13H4.3fdundefinedFtank1.stl","affiliations":{"Deck":"WheelHouseTop","SFI":"101"},"referenceState":{"xCentre":2,"yCentre":0,"zBase":4}},{"id":"Tank2","baseObject":"TankL2.1B6.13H4.3fdundefinedFundefined","affiliations":{"Deck":"BridgeDeck","SFI":"102"},"referenceState":{"xCentre":5,"yCentre":0,"zBase":2}},{"id":"Tank3","baseObject":"TankL4.199999999999999B6.7H4.3fdundefinedFtank2.stl","affiliations":{"Deck":"BridgeDeck","SFI":"103"},"referenceState":{"xCentre":14,"yCentre":0,"zBase":2}}]
+
+}
\ No newline at end of file
diff --git a/source/CoreClasses/Ship.js b/source/CoreClasses/Ship.js
index 1885542..0414559 100644
--- a/source/CoreClasses/Ship.js
+++ b/source/CoreClasses/Ship.js
@@ -83,6 +83,7 @@ Object.assign(Ship.prototype, {
 		console.info("Calculated draft: %.2f", t);
 		return t;
 	},
+	//Should separate between longitudinal and transverse GM too
     calculateStability(shipState){
         let T = this.calculateDraft(shipState);
         let ha = this.structure.hull.calculateAttributesAtDraft(T);
@@ -94,7 +95,7 @@ Object.assign(Ship.prototype, {
             vol = Lwl * B * T * cb;
         }
         let KG = this.getWeight(shipState).cg.z;
-        let I = ha.Iywp * 1000;
+        let I = ha.Iywp;
         let KB = 0.52 * T;
         let BM = I / vol;
         let GM = KB + BM - KG;

From 2eec29e45a4fd61f1ccaf2a8b2c227e6a276039d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Vicente=20Alejandro=20Iv=C3=A1=C3=B1ez=20Encinas?=
 
Date: Mon, 23 Oct 2017 18:21:24 +0200
Subject: [PATCH 03/20] Add some values and correct existing ones

---
 examples/blockCase_compare.html               | 103 +++++++++++++-----
 .../data/ship_specifications/blockCase.json   |   2 +-
 2 files changed, 75 insertions(+), 30 deletions(-)

diff --git a/examples/blockCase_compare.html b/examples/blockCase_compare.html
index 00315c3..b610827 100644
--- a/examples/blockCase_compare.html
+++ b/examples/blockCase_compare.html
@@ -122,6 +122,24 @@ 
Comparation

                     
                     
                     
+                        
+                           Draft Amidships (m) 
+                            
+                            
+                            
+                        
+                        
+                           Draft at FP (m) 
+                            
+                            
+                            
+                        
+                        
+                           Draft AP (m) 
+                            
+                            
+                            
+                        
                         
                            Displacement (t) 
                             
@@ -131,14 +149,14 @@ 
Comparation

                         
                            Waterline Lenght (m) 
                             
-                            
-                            
+                            
+                            
                         
                         
                            Maximum Waterline Breath (m) 
                             
-                            
-                            
+                            
+                            
                         
                         
                            Wetted area (m2) 
@@ -149,8 +167,8 @@ 
Comparation

                         
                            Waterplane Area (m2) 
                             
-                            
-                            
+                            
+                            
                         
                         
                            Cb 
@@ -289,50 +307,77 @@ 
Comparation

             document.getElementById("GM2").innerHTML = GMc;
             
             //hand values
-            let disp = 111;
+            let disp = 102.4;
+            let draftms = 1.989;
+            let draftap = 1.98;
+            let draftfp = 1.999;
             let lwl = 20;
-            let bwl = 5.292;
-            let Awet = 144.091;
-            let Awp = 104.059;
-            let Cb = 0.975;
-            let LCBx = 9.885;
-            let LCFx = 9.943;
-            let KB = 1.388;
-            let KG = 3.473;
-            let BM = 2.168;
-            let GM = 0.083;
+            let bwl = 5.046;
+            let Awet = 137.989;
+            let Awp = 99.95;
+            let Cb = disp / (draftms * lwl * bwl * 1.025);
+            let LCBx = 10.065;
+            let LCFx = 10.032;
+            let KB = 1.333;
+            let KG = 3.513;
+            let BM = 2.082;
+            let GM = KB + BM - KG;
+            let trim = -0.019;
             
-            document.getElementById("disp").innerHTML = disp;
-            document.getElementById("lwl").innerHTML = lwl;
-            document.getElementById("bwl").innerHTML = bwl;
-            document.getElementById("Awet").innerHTML = Awet;
+            document.getElementById("disp").innerHTML = disp.toFixed(3);
+            document.getElementById("draftms").innerHTML = draftms.toFixed(3);
+            document.getElementById("draftfp").innerHTML = draftfp.toFixed(3);
+            document.getElementById("draftap").innerHTML = draftap.toFixed(3);
+            document.getElementById("lwl").innerHTML = lwl.toFixed(3);
+            document.getElementById("bwl").innerHTML = bwl.toFixed(3);
+            document.getElementById("Awet").innerHTML = Awet.toFixed(3);
             document.getElementById("Awp").innerHTML = Awp;
-            document.getElementById("Cb").innerHTML = Cb;
-            document.getElementById("LCBx").innerHTML = LCBx;
-            document.getElementById("LCFx").innerHTML = LCFx;
-            document.getElementById("KB").innerHTML = KB;
-            document.getElementById("KG").innerHTML = KG;
-            document.getElementById("BM").innerHTML = BM;
-            document.getElementById("GM").innerHTML = GM;
+            document.getElementById("Cb").innerHTML = Cb.toFixed(3);
+            document.getElementById("LCBx").innerHTML = LCBx.toFixed(3);
+            document.getElementById("LCFx").innerHTML = LCFx.toFixed(3);
+            document.getElementById("KB").innerHTML = KB.toFixed(3);
+            document.getElementById("KG").innerHTML = KG.toFixed(3);
+            document.getElementById("BM").innerHTML = BM.toFixed(3);
+            document.getElementById("GM").innerHTML = GM.toFixed(3);
             
             //calculated values
             let KBc = a.KB.toFixed(3);
             let KGc = a.KG.toFixed(3);
             let BMc = a.BM.toFixed(3);
+            let draftc = ship.calculateDraft().toFixed(3);
+            let lwlc = ship.structure.hull.calculateAttributesAtDraft().LWL;
+            let bwlc = ship.structure.hull.calculateAttributesAtDraft().BWL;
+            let Awpc = ship.structure.hull.calculateAttributesAtDraft().Awp;
+            
             document.getElementById("KBc").innerHTML = KBc;
             document.getElementById("KGc").innerHTML = KGc;
             document.getElementById("BMc").innerHTML = BMc;
+            document.getElementById("draftc").innerHTML = draftc;
+            document.getElementById("lwlc").innerHTML = lwlc.toFixed(3);
+            document.getElementById("bwlc").innerHTML = bwlc.toFixed(3);
+            document.getElementById("Awpc").innerHTML = Awpc.toFixed(3);
             
             //% diff
             let KBd = (1-KB/KBc)*100;
             let KGd = (1-KG/KGc)*100;
             let BMd = (1-BM/BMc)*100;
             let GMd = (1-GM/GMc)*100;
+            let draftd = (1-draftms/draftc)*100;
+            let lwld = (1-lwl/lwlc)*100;
+            let bwld = (1-bwl/bwlc)*100;
+            let awpd = (1-Awp/Awpc)*100;
+            
             document.getElementById("KBd").innerHTML = KBd.toFixed(3);
             document.getElementById("KGd").innerHTML = KGd.toFixed(3);
             document.getElementById("BMd").innerHTML = BMd.toFixed(3);
             document.getElementById("GMd").innerHTML = GMd.toFixed(3);
-
+            document.getElementById("draftd").innerHTML = draftd.toFixed(3);
+            document.getElementById("lwld").innerHTML = lwld.toFixed(3);
+            document.getElementById("bwld").innerHTML = bwld.toFixed(3);
+            document.getElementById("awpd").innerHTML = awpd.toFixed(3);
+            
+            
+            console.log(ship.structure.hull.calculateAttributesAtDraft())
 		}
 		
 		function animate() {
diff --git a/examples/data/ship_specifications/blockCase.json b/examples/data/ship_specifications/blockCase.json
index b3a4601..ee5e741 100644
--- a/examples/data/ship_specifications/blockCase.json
+++ b/examples/data/ship_specifications/blockCase.json
@@ -5,7 +5,7 @@
 	"calculationParameters": {
 		"LWL_design": 20,
 		"Draft_design": 2,
-		"Cb_design": 1,
+		"Cb_design": 0.5,
 		"speed": 12,
 	    "crew" : 20,
 		"K" : 0.032,

From 2a0610fdef462da5099007c17f7383a2a2fe8816 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Vicente=20Alejandro=20Iv=C3=A1=C3=B1ez=20Encinas?=
 
Date: Mon, 23 Oct 2017 18:26:30 +0200
Subject: [PATCH 04/20] Add Cb

---
 examples/blockCase_compare.html | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/examples/blockCase_compare.html b/examples/blockCase_compare.html
index b610827..720e381 100644
--- a/examples/blockCase_compare.html
+++ b/examples/blockCase_compare.html
@@ -173,7 +173,7 @@ 
Comparation

                         
                            Cb 
                             
-                            
+                            
                             
                         
                         
@@ -348,6 +348,7 @@ 
Comparation

             let lwlc = ship.structure.hull.calculateAttributesAtDraft().LWL;
             let bwlc = ship.structure.hull.calculateAttributesAtDraft().BWL;
             let Awpc = ship.structure.hull.calculateAttributesAtDraft().Awp;
+            let cbc = ship.structure.hull.calculateAttributesAtDraft().Cb;
             
             document.getElementById("KBc").innerHTML = KBc;
             document.getElementById("KGc").innerHTML = KGc;
@@ -356,6 +357,7 @@ 
Comparation

             document.getElementById("lwlc").innerHTML = lwlc.toFixed(3);
             document.getElementById("bwlc").innerHTML = bwlc.toFixed(3);
             document.getElementById("Awpc").innerHTML = Awpc.toFixed(3);
+            document.getElementById("cbc").innerHTML = cbc.toFixed(3);
             
             //% diff
             let KBd = (1-KB/KBc)*100;

From d2b18d44f7a4f0cb0b9172303b64d146d91944bf Mon Sep 17 00:00:00 2001
From: EliasHasle 
Date: Tue, 24 Oct 2017 22:57:27 +0200
Subject: [PATCH 05/20] Added check for missing draft argument to
 Hull.prototype.calculateAttributesAtDraft

Also improved comments to highlight that an argument is required.
---
 build/Vessel.js            | 9 ++++++++-
 source/CoreClasses/Hull.js | 7 +++++++
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/build/Vessel.js b/build/Vessel.js
index c684278..75e7bce 100644
--- a/build/Vessel.js
+++ b/build/Vessel.js
@@ -1,4 +1,4 @@
-//Vessel.js library, built 2017-10-18 09:13:35.951607, Checksum: 98f5cf6e023f4efb329738cead04ff77
+//Vessel.js library, built 2017-10-24 22:53:41.412445, Checksum: c9f7fffc8d35ee190cfc25d553418214
 /*
 Import like this in HTML:
 
@@ -1116,6 +1116,7 @@ Object.assign(Hull.prototype, {
 	},
 	//Unoptimized, some redundant repetitions of calculations.
 	//NOT DONE YET. Outputs lots of NaN values.
+	//Important: calculateAttributesAtDraft takes one mandatory parameter T. (The function defined here is immediately called during construction of the prototype, and returns the proper function.)
 	calculateAttributesAtDraft: function() {
 		function levelCalculation(hull,
 			z,
@@ -1188,7 +1189,13 @@ Object.assign(Hull.prototype, {
 			return lev;
 		}
 		
+		//Here is the returned function calculateAttributesAtDraft(T):
 		return function(T) {
+			if (T === undefined) {
+				console.error("Hull.prototype.calculateAttributesAtDraft(T): No draft specified. Returning undefined.");
+				return;
+			}
+			
 			let wls = this.halfBreadths.waterlines.map(wl=>this.attributes.Depth*wl);
 			
 			//This is the part that can be reused as long as the geometry remains unchanged:
diff --git a/source/CoreClasses/Hull.js b/source/CoreClasses/Hull.js
index d7c585d..a8fa586 100644
--- a/source/CoreClasses/Hull.js
+++ b/source/CoreClasses/Hull.js
@@ -309,6 +309,7 @@ Object.assign(Hull.prototype, {
 	},
 	//Unoptimized, some redundant repetitions of calculations.
 	//NOT DONE YET. Outputs lots of NaN values.
+	//Important: calculateAttributesAtDraft takes one mandatory parameter T. (The function defined here is immediately called during construction of the prototype, and returns the proper function.)
 	calculateAttributesAtDraft: function() {
 		function levelCalculation(hull,
 			z,
@@ -381,7 +382,13 @@ Object.assign(Hull.prototype, {
 			return lev;
 		}
 		
+		//Here is the returned function calculateAttributesAtDraft(T):
 		return function(T) {
+			if (T === undefined) {
+				console.error("Hull.prototype.calculateAttributesAtDraft(T): No draft specified. Returning undefined.");
+				return;
+			}
+			
 			let wls = this.halfBreadths.waterlines.map(wl=>this.attributes.Depth*wl);
 			
 			//This is the part that can be reused as long as the geometry remains unchanged:

From 21d76b6381e6fe269558f94ddb896f3d9a07f4b4 Mon Sep 17 00:00:00 2001
From: EliasHasle 
Date: Tue, 24 Oct 2017 23:22:22 +0200
Subject: [PATCH 06/20] Sketched comparison with rendered json instead of
 handmade tables

---
 examples/blockCase_compare_json.html | 86 ++++++++++++++++++++++++++++
 1 file changed, 86 insertions(+)
 create mode 100644 examples/blockCase_compare_json.html

diff --git a/examples/blockCase_compare_json.html b/examples/blockCase_compare_json.html
new file mode 100644
index 0000000..92110e2
--- /dev/null
+++ b/examples/blockCase_compare_json.html
@@ -0,0 +1,86 @@
+
+
+
+
+	Comparison of manually calculations and library calculations
+	
+	
+	
+
+
+
    
+
+
+
+
\ No newline at end of file

From 132a3dbc44590029e17c3d7cd0b751b810d48770 Mon Sep 17 00:00:00 2001
From: EliasHasle 
Date: Wed, 25 Oct 2017 14:16:48 +0200
Subject: [PATCH 07/20] Added links to more test applications

---
 examples/index.html | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/examples/index.html b/examples/index.html
index c44ceb5..e977d40 100644
--- a/examples/index.html
+++ b/examples/index.html
@@ -27,7 +27,11 @@ 
    Examples
    
 
    
 
    Ship3D_with_pretty_JSON.html
     
    A 3D visualization of the provided ship specification, along with the specification data in a format with expand and collapse.
    
 
    Ship3D_test.html
     
    A fullscreen 3D visualization of the provided ship specification.
    
+
    Ship3D_load_state_controls.html
     
    Very simple test of changing tank fullness with a slider to see change of draft.
    
 
    Vessel.js_test.html
     
    A script that does a few calculations based on the provided specification, and outputs to the screen and the console.
    
+
    blockCase_compare.html
     
    Comparison between manual and automatic calculations on a very simple ship.
    
+
    blockCase_compare_json.html
     
    Sketched alternative implementation of same comparison.
    
+
    PX121_calculations.html
     
    Some calculations on a test specification partially based on an Ulstein vessel.
    
 
    
 
 

From c095abb5539587d09eb3fb2c6753539191690027 Mon Sep 17 00:00:00 2001
From: EliasHasle 
Date: Wed, 25 Oct 2017 16:48:49 +0200
Subject: [PATCH 08/20] Minor improvements to comparison (JSON version)

---
 examples/blockCase_compare_json.html | 58 ++++++++++++++++++++++------
 1 file changed, 47 insertions(+), 11 deletions(-)

diff --git a/examples/blockCase_compare_json.html b/examples/blockCase_compare_json.html
index 92110e2..ebf834f 100644
--- a/examples/blockCase_compare_json.html
+++ b/examples/blockCase_compare_json.html
@@ -27,10 +27,12 @@
 function compToScreen(s, man, machine) {
 	let log = document.getElementById("log");
 	let entry = document.createElement("li");
-	entry.innerHTML = s;
+	entry.innerHTML = s+" ";
 
 	let left = document.createElement("div");
 	left.style.display = "inline-block";
+	left.style.verticalAlign = "top";
+	left.appendChild(document.createTextNode("Man"));
 	let preLeft = document.createElement("pre");
 	preLeft.appendChild(renderjson(man));
 	left.appendChild(preLeft);
@@ -38,6 +40,8 @@
 	
 	let right = document.createElement("div");
 	right.style.display = "inline-block";
+	right.style.verticalAlign = "top";
+	right.appendChild(document.createTextNode("Machine"));
 	let preRight = document.createElement("pre");
 	preRight.appendChild(renderjson(machine));
 	right.appendChild(preRight);
@@ -54,26 +58,58 @@
 function useShip(ship) {
 	v = ship; //just to have an entry point for debugging.
 
-	renderjson.set_show_to_level(1);
+	renderjson.set_show_to_level(2);
+	renderjson.set_sort_objects(true);
 	
 	logToScreen("Ship: ", ship);
 	//logToScreen("Design state: ", ship.designState);
 
 	//Hand-calculated values from @MrEranwe
+	//Some of these values are not calculated by the library yet.
+	//And the library does not account for trim when calculating Cb etc.
+	let disp = 102.4;
+	let draftms = 1.989;
+	//let draftap = 1.98;
+	//let draftfp = 1.999;
+	let lwl = 20;
+	let bwl = 5.046;
+	let Awet = 137.989;
+	let Awp = 99.95;
+	let Cb = disp / (draftms * lwl * bwl * 1.025);
+	let LCBx = 10.065;
+	let LCFx = 10.032;
+	let KB = 1.333;
+	let KG = 3.513;
+	let BM = 2.082;
+	let GM = KB + BM - KG;
+	//let trim = -0.019;
 	
 	let w = ship.getWeight();
-	compToScreen("Weight at design state:", {}, w);
+	compToScreen("Weight at design state:", {
+		mass: 1000*disp,
+		cg: {
+			x: undefined,
+			y: undefined,
+			z: KG
+		}
+	}, w);
 	
 	let T = ship.calculateDraft();
-	logToScreen("Calculated draft: Man: 1.989, Machine: " + T);
-	
-	let ha = v.structure.hull.calculateAttributesAtDraft(T);
+	logToScreen("Calculated draft: Man: "+ draftms+ ", Machine: "+ T);
+	logToScreen("Assuming draft = " + draftms + ":");
+	let ha = v.structure.hull.calculateAttributesAtDraft(draftms);
 	compToScreen("Calculated hull values at calculated draft: ", {
-		LWL: 20,
-		BWL: 5.046,
-		As: 137.989,
-		Awp: 99.95,
-		Cb: 0.498
+		LWL: lwl,
+		BWL: bwl,
+		As: Awet,
+		Awp: Awp,
+		Cb: Cb,
+		Cv: {
+			x: LCBx,
+			y: undefined,
+			z: KB
+		},
+		xcwp: LCFx
 		//more
 	}, ha);
 	

From 2b24ae87c667a8f0796d65c59c15f9865472bd78 Mon Sep 17 00:00:00 2001
From: EliasHasle 
Date: Thu, 26 Oct 2017 13:35:03 +0200
Subject: [PATCH 09/20] Made blockCase compare more useful

Now the differences are easy to spot, even though it does not highlight
errors. The new manual calculations ignore trim. What is currently
missing here, is weight calculation and stability calculations. But for
now it should be enough to focus on getting the geometric calculations
right.
---
 examples/blockCase_compare_json.html | 83 ++++++++++++++++++++++------
 1 file changed, 66 insertions(+), 17 deletions(-)

diff --git a/examples/blockCase_compare_json.html b/examples/blockCase_compare_json.html
index ebf834f..a21cf53 100644
--- a/examples/blockCase_compare_json.html
+++ b/examples/blockCase_compare_json.html
@@ -58,16 +58,17 @@
 function useShip(ship) {
 	v = ship; //just to have an entry point for debugging.
 
-	renderjson.set_show_to_level(2);
 	renderjson.set_sort_objects(true);
-	
+		
 	logToScreen("Ship: ", ship);
 	//logToScreen("Design state: ", ship.designState);
 
+	renderjson.set_show_to_level("all");
+		
 	//Hand-calculated values from @MrEranwe
 	//Some of these values are not calculated by the library yet.
 	//And the library does not account for trim when calculating Cb etc.
-	let disp = 102.4;
+	/*let disp = 102.4;
 	let draftms = 1.989;
 	//let draftap = 1.98;
 	//let draftfp = 1.999;
@@ -82,36 +83,84 @@
 	let KG = 3.513;
 	let BM = 2.082;
 	let GM = KB + BM - KG;
-	//let trim = -0.019;
+	//let trim = -0.019;*/
 	
 	let w = ship.getWeight();
-	compToScreen("Weight at design state:", {
+	logToScreen("Assuming weight from library calculation: ", w);
+	//Weight cannot be compared until the method for calculating it has been specified. Right now the library method is faulty anyway.
+	/*compToScreen("Weight at design state: ", {
 		mass: 1000*disp,
 		cg: {
 			x: undefined,
 			y: undefined,
 			z: KG
 		}
-	}, w);
+	}, w);*/
+
+	
+	//Simplified calculations for triangular prism case, without trim:
+	//Assuming that T=Depth) {
+		logToScreen("Invalid calculations, because T>=Depth. T=" + T);
+	}
+	let BWL = (BOA/Depth)*T;
+	let LWL = LOA;
+	let Vs = 0.5*LWL*BWL*T;
+	//let Vbb = BWL*LWL*T;// = LOA*(BOA/Depth)*T^2
+	let Cb = 0.5;// = Vs/Vbb
+	let Ap = 0.5*BWL*T;
+	let As = 2*Ap + 2*LWL*Math.sqrt(T**2 + (0.5*BWL)**2);
+	let Awp = LWL*BWL;
+	let Cwp = 1;
+	let Cv = {x:0.5*LOA, y:0, z:(2/3.)*T}
+	let Ixwp = LWL*BWL**3/12.;
+	let Iywp = BWL*LWL**3/12.;
+	let xcwp = 0.5*LWL;
+	let ycwp = 0;
+	let minXs = 0;
+	let maxXs = LWL;
+	let minYs = -0.5*BWL;
+	let maxYs = 0.5*BWL;
+	let LBP = LWL-ha.APP;
 	
-	let T = ship.calculateDraft();
-	logToScreen("Calculated draft: Man: "+ draftms+ ", Machine: "+ T);
-	logToScreen("Assuming draft = " + draftms + ":");
-	let ha = v.structure.hull.calculateAttributesAtDraft(draftms);
+	let Tc = ship.calculateDraft();
+	logToScreen("Calculated draft: Man: "+ T + ", Machine: "+ Tc);
+	logToScreen("Assuming draft = " + T + ":");
+	let hca = v.structure.hull.calculateAttributesAtDraft(T);
 	compToScreen("Calculated hull values at calculated draft: ", {
-		LWL: lwl,
-		BWL: bwl,
-		As: Awet,
+		LWL: LWL,
+		BWL: BWL,
+		Ap: Ap,
+		As: As,
 		Awp: Awp,
+		Cwp: Cwp,
 		Cb: Cb,
-		Cv: {
+		Cv: Cv,/*{
 			x: LCBx,
 			y: undefined,
 			z: KB
-		},
-		xcwp: LCFx
+		},*/
+		xcwp: xcwp,//LCFx
+		ycwp: ycwp,
+		Ixwp: Ixwp,
+		Iywp: Iywp,
+		Vs: Vs,
+		//Vbb: Vbb,
+		minXs: minXs,
+		maxXs: maxXs,
+		minYs: minYs,
+		maxYs: maxYs,
+		LBP: LBP
 		//more
-	}, ha);
+	}, hca);
 	
 	//let wl = v.structure.hull.getWaterline(T);
 	//logToScreen("Waterline at calculated draft: ", wl);

From a648ef46dc2860b24f0c419b93fe28bd5b9a37d0 Mon Sep 17 00:00:00 2001
From: EliasHasle 
Date: Thu, 26 Oct 2017 22:25:50 +0200
Subject: [PATCH 10/20] Correctd Ixwp and Iywp, and some more things

---
 build/Vessel.js                              |  161 +-
 examples/blockCase_compare_json.html         |    3 +-
 examples/data/ship_specifications/PX121.json | 4112 +++++++++++++++++-
 source/CoreClasses/Hull.js                   |   93 +-
 source/CoreClasses/Ship.js                   |    2 +-
 source/CoreClasses/Structure.js              |    6 +-
 source/functions/areaCalculations.js         |   14 +-
 source/functions/interpolation.js            |   32 +-
 source/functions/volumeCalculations.js       |   12 +-
 9 files changed, 4221 insertions(+), 214 deletions(-)

diff --git a/build/Vessel.js b/build/Vessel.js
index 75e7bce..a9efb52 100644
--- a/build/Vessel.js
+++ b/build/Vessel.js
@@ -1,4 +1,4 @@
-//Vessel.js library, built 2017-10-24 22:53:41.412445, Checksum: c9f7fffc8d35ee190cfc25d553418214
+//Vessel.js library, built 2017-10-26 22:24:21.268970, Checksum: 69efcb4dadd22b34bfc962e870f528a7
 /*
 Import like this in HTML:
 
@@ -63,22 +63,6 @@ function crossProduct(u,v) {
 
 //Some interpolation helpers. Only linear and bilinear for now.
 
-//linear interpolation
-//Defaults are not finally decided
-//returns NaN if a and b are NaN or mu is NaN.
-function lerp(a, b, mu=0.5) {
-	if (isNaN(a)) return b;
-	if (isNaN(b)) return a;
-	return (1-mu)*a+mu*b;
-}
-
-//Test. Not safe yet.
-function linearFromArrays(xx, yy, x) {
-	let {index, mu} = bisectionSearch(xx, x);
-	if (index === undefined || mu === undefined) return 0;
-	return lerp(yy[index], yy[index+1], mu);	
-}
-
 /*Function that takes a sorted array as input, and finds the last index that holds a numerical value less than, or equal to, a given value.
 Returns an object with the index and an interpolation parameter mu that gives the position of value between index and index+1.
 */
@@ -105,6 +89,22 @@ function bisectionSearch(array, value) {
 	return {index, mu};
 }
 
+//linear interpolation
+//Defaults are not finally decided
+//returns NaN if a and b are NaN or mu is NaN.
+function lerp(a, b, mu=0.5) {
+	if (isNaN(a)) return b;
+	if (isNaN(b)) return a;
+	return (1-mu)*a+mu*b;
+}
+
+//Test. Not safe yet.
+function linearFromArrays(xx, yy, x) {
+	let {index, mu} = bisectionSearch(xx, x);
+	if (index === undefined || mu === undefined) return 0;
+	return lerp(yy[index], yy[index+1], mu);	
+}
+
 function bilinearUnitSquareCoeffs(z00, z01, z10, z11) {
 	let a00 = z00;
 	let a10 = z10-z00;
@@ -159,7 +159,7 @@ function bilinear(x1, x2, y1, y2, z11, z12, z21, z22, x, y) {
 
 //All inputs are numbers. The axes are given by a single coordinate.
 function steiner(I, A, sourceAxis, targetAxis) {
-	return I + A*(sourceAxis-targetAxis)^2;
+	return I + A*(sourceAxis-targetAxis)**2;
 }
 
 //Calculate area, center, Ix, Iy.
@@ -174,18 +174,18 @@ function trapezoidCalculation(xbase0, xbase1, xtop0, xtop1, ybase, ytop) {
 	let yc = (a==0 && b==0) ? ybase+0.5*h : ybase + h*(2*a+b)/(3*(a+b));
 	let d = xbase0+0.5*a; //shorthand
 	let xc = h===0 ? 0.25*(xbase0+xbase1+xtop0+xtop1) : d + (xtop0+0.5*b-d)*(yc-ybase)/h;
-	let Ix = (a==0 && b== 0) ? 0 : h^3*(a^2+4*a*b+b^2)/(36*(a+b));
+	let Ix = (a==0 && b== 0) ? 0 : h**3*(a**2+4*a*b+b**2)/(36*(a+b));
 
 	//For Iy I must decompose (I think negative results will work fine):
 	let Art1 = 0.5*(xtop0-xbase0)*h;
 	let xcrt1 = xbase0 + (xtop0-xbase0)/3;
-	let Iyrt1 = (xtop0-xbase0)^3*h/36;
+	let Iyrt1 = (xtop0-xbase0)**3*h/36;
 	let Arec = (xbase1-xtop0)*h;
 	let xcrec = 0.5*(xtop0+xbase1);
-	let Iyrec = (xbase1-xtop0)^3*h/12;
+	let Iyrec = (xbase1-xtop0)**3*h/12;
 	let Art2 = 0.5*(xbase1-xtop1)*h;
 	let xcrt2 = (xtop1 + (xbase1-xtop1)/3);
-	let Iyrt2 = (xbase1-xtop1)^3*h/36;
+	let Iyrt2 = (xbase1-xtop1)**3*h/36;
 
 	let Iy = steiner(Iyrt1, Art1, xcrt1, xc)
 		+ steiner(Iyrec, Arec, xcrec, xc)
@@ -242,10 +242,8 @@ function combineAreas(array) {
 
 //x and y here refers to coordinates in the plane that is being calculated on.
 function sectionCalculation({xs, ymins, ymaxs}) {
-	console.groupCollapsed("sectionCalculation");
+	console.group/*Collapsed*/("sectionCalculation");
 	console.info("Arguments (xs, ymins, ymaxs): ", arguments[0]);
-	
-	//Needed for Cwp (not a very efficient calculation, maybe):
 
 	let calculations = [];
 	for (let i = 0; i < xs.length-1; i++) {
@@ -278,19 +276,19 @@ function bilinearPatchColumnCalculation(x1, x2, y1, y2, z11, z12, z21, z22) {
 	= INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x) dx]
 	= a00 + 0.5*a10 + 0.5*a01 + 0.25*a11
 	*/
-	let Ab = X*Y;
+	let Ab = X*Y; //area of base of patch column
 	let zAvg = (a00 + 0.5*a10 + 0.5*a01 + 0.25*a11);
 	let V = Math.abs(Ab*zAvg); //new: absolute value
 	let zc = 0.5*zAvg;
 	/*
-	To find xc, I need to integrate x*z over the unit square, and scale and translate to world coordinates afterwards:
+	To find xc, I need to integrate x*z over the unit square, and scale and translate to ship coordinates afterwards:
 	INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x)*x dx]
 	= 0.5*a00 + a10/3 + 0.25*a01 + a11/6
 	Scale and translate:*/
-	let xc = y1 + X*(0.5*a00 + a10/3 + 0.25*a01 + a11/6)
+	let xc = x1 + X*(0.5*a00 + a10/3 + 0.25*a01 + a11/6);
 	
 	//Similar for yc:
-	let yc = y1 + Y*(0.5*a00 + 0.25*a10 + a01/3 + a11/6)
+	let yc = y1 + Y*(0.5*a00 + 0.25*a10 + a01/3 + a11/6);
 	
 	//new: absolute value (OK?)
 	let As = Math.abs(bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22));
@@ -330,6 +328,8 @@ function combineVolumes(array) {
 }
 
 //For wetted area. I think this is right, but it is not tested.
+//The numerical integral will not scale well with larger geometries.
+//Then the full analytical solution is needed.
 function bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22, segs=10) {
 	let [b00,b10,b01,b11] = bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22);
 	/*
@@ -343,6 +343,8 @@ function bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22, segs=10) {
 	Tx X Ty = (-(b10+b11*y), -(b01+b11*x), 1)
 	|Tx X Ty| = sqrt((b10+b11*y)^2 + (b01+b11*x)^2 + 1)
 	
+	Now, to get the area I need to integrate |Tx X Ty| over X,Y.
+	
 	Wolfram Alpha gave me this for the inner integral using x (indefinite):
 	integral sqrt((b01 + b11 x)^2 + 1 + (b10+b11*y)^2) dx = ((b01 + b11*x) sqrt((b01 + b11*x)^2 + 1 + (b10+b11*y)^2) + (1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x))/(2*b11) + constant
 	That means this for the definite integral:
@@ -624,7 +626,7 @@ Object.assign(Ship.prototype, {
 	constructor: Ship,
 	setFromSpecification: function(specification) {
 		this.attributes = specification.attributes || {};
-		this.structure = new Structure(specification.structure,this);
+		this.structure = new Structure(specification.structure/*,this*/);
 		//baseObjects and derivedObjects are arrays in the specification, but are objects (hashmaps) in the constructed ship object:
 		this.baseObjects = {};
 		for (let i = 0; i < specification.baseObjects.length; i++) {
@@ -710,8 +712,8 @@ Object.assign(Ship.prototype, {
 	}
 });//@EliasHasle
 
-function Structure(spec, ship) {
-	this.ship = ship;
+function Structure(spec/*, ship*/) {
+	//this.ship = ship;
 	JSONSpecObject.call(this, spec);
 }
 Structure.prototype = Object.create(JSONSpecObject.prototype);
@@ -773,7 +775,7 @@ Object.assign(Structure.prototype, {
 			let zc = d.zFloor+0.5*d.thickness;
 			let yc = d.yCentre;
 			let b = d.breadth;
-			let wlc = this.hull.waterlineCalculation(zc, {minX: d.xAft, maxX: d.xFwd, minY: yc-0.5*b, maxY: yc+0.5*b});
+			let wlc = this.hull.waterlineCalculation(zc, {minX: d.xAft, maxX: d.xFwd, minY: yc-0.5*b, maxY: yc+0.5*b}, 3);
 			components.push({
 				//approximation
 				mass: wlc.Awp*d.thickness*d.density,
@@ -856,9 +858,11 @@ Object.assign(Hull.prototype, {
 		return output;
 	},
 	/*
+	This method is a mess, and introducing three different NaN correction has made it even worse in some ways. None of the correction modes fit for the usual cases. Maybe we should have a default mode that extrapolates to top NaNs and sets all other NaNs to zero. UPDATE: Now this is implemented as mode 3. Maybe we can skip the other ones.
+	
 	Input:
 	z: level from bottom of ship (absolute value in meters)
-	nanCorrectionMode: 0 to set all NaNs to zero, 1 to output NaNs, set all NaNs to zero, 2 to replace NaNs with interpolated or extrapolated values.	
+	nanCorrectionMode: 0 to set all NaNs to zero, 1 to output NaNs, 2 to replace NaNs with interpolated or extrapolated values, 3 to extrapolate to top NaNs and set all other NaNs to zero.
 	*/
 	getWaterline: function(z, nanCorrectionMode=1) {
 		let ha = this.attributes;
@@ -870,7 +874,7 @@ Object.assign(Hull.prototype, {
 		let {index: a, mu: mu} = bisectionSearch(wls, zr);
 		let wl;
 		if (a<0) {
-			if (nanCorrectionMode===0) {
+			if (nanCorrectionMode===0 || nanCorrectionMode===3) {
 				console.warn("getWaterLine: z below lowest defined waterline. Defaulting to zeros.");
 				return new Array(sts.length).fill(0);
 			}
@@ -884,7 +888,7 @@ Object.assign(Hull.prototype, {
 				mu=0;
 				//wl = tab[a].slice();
 			}			
-		} else if (a>/*=*/wls.length-1) {
+		} else if (a >= wls.length-1) {
 			if (nanCorrectionMode===0) {
 				console.warn("getWaterLine: z above highest defined waterline. Defaulting to zeros.");
 				return new Array(sts.length).fill(0);
@@ -893,7 +897,7 @@ Object.assign(Hull.prototype, {
 				console.warn("getWaterLine: z above highest defined waterline. Outputting NaNs.");
 				return new Array(sts.length).fill(null);
 			}
-			else /*nanCorrectionMode===2*/ {
+			else /*nanCorrectionMode===2 || nanCorrectionMode===3*/ {
 				console.warn("getWaterLine: z above highest defined waterline. Proceeding with highest data entry.");
 				a = wls.length-2; //if this level is defined...
 				mu=1;
@@ -906,7 +910,7 @@ Object.assign(Hull.prototype, {
 		for (let j = 0; j < wl.length; j++) {
 			if (nanCorrectionMode === 0) {
 				if (a+1 > wls.length-1) {
-					wl[j] = lerp(tab[a][j], 0, 0.5);
+					wl[j] = lerp(tab[a][j], 0, 0.5); //Suspicious!
 				} else {
 					wl[j] = lerp(tab[a][j] || 0, tab[a+1][j] || 0, mu || 0.5);
 				}
@@ -916,7 +920,7 @@ Object.assign(Hull.prototype, {
 				} else {
 					wl[j] = lerp(tab[a][j], tab[a+1][j], mu);
 				}
-			} else {
+			} else if (nanCorrectionMode === 2) {
 				//If necessary, sample from below
 				let b = a;
 				while (b>0 && isNaN(tab[b][j])) {
@@ -948,14 +952,62 @@ Object.assign(Hull.prototype, {
 						upper = tab[c][j];
 					}
 				}
-				mu = c===b ? 0 : (a+(mu||0.5)-b)/(c-b);
+				mu = c===b ? 0 : (a+(mu||0.5)-b)/(c-b); //what is this? Some kind of self-explanatory code, perhaps?
 				wl[j] = lerp(lower, upper, mu);
+			} else /*nanCorrectionMode === 3*/ {
+				let lower, upper;
+				let b;
+				if (!isNaN(tab[a][j])) {
+					lower = tab[a][j];					
+				} else {
+					b = a+1;
+					while(b < wls.length && isNaN(tab[b][j])) {
+						b++;
+					}
+					if (b !== wls.length) {
+						//Inner NaN
+						lower = 0;
+					} else {
+						//Upper NaN, search below:
+						b = a-1;
+						while (b >= 0 && isNaN(tab[b][j])) {
+							b--;
+						}
+						if (b===-1) {
+							//No number found:
+							lower = 0;
+							upper = 0;
+						} else {
+							lower = tab[b][j];
+							upper = lower;
+						}
+					}
+				}
+				if (upper !== undefined) {}
+				else if (!isNaN(tab[a+1][j])) {
+					upper = tab[a+1][j];
+				} else {
+					//The cell value is NaN.
+					//Upper is not defined.
+					//That means either tab[a][j] is a number
+					//or tab[a][j] is an inner NaN and
+					//there exists at least one number above it.
+					//In both cases I have to check above a+1.
+					b = a+2;
+					while (b < wls.length && isNaN(tab[b][j])) {
+						b++;
+					}
+					if (b === wls.length) upper = lower;
+					else {
+						upper = tab[b][j];
+					}
+				}
+				wl[j] = lerp(lower, upper, mu || 0.5);
 			}
 			
 			//Scale numerical values
 			if (!isNaN(wl[j])) wl[j] *= 0.5*ha.BOA;
 		}
-
 		return wl;
 	},
 	getStation: function(x) {
@@ -991,14 +1043,14 @@ Object.assign(Hull.prototype, {
 	
 	//THIS is a candidate for causing wrong Ix, Iy values.
 	//Much logic that can go wrong.
-										//typically deck bounds
-	waterlineCalculation: function(z, bounds) {
+									//typically deck bounds
+	waterlineCalculation: function(z, bounds, nanCorrectionMode=3) {
 		let {minX, maxX, minY, maxY} = bounds || {};
 
-		console.groupCollapsed("waterlineCalculation.");
-		console.info("Arguments: z=", z, " Boundaries: ", arguments[1]);
+		console.group/*Collapsed*/("waterlineCalculation.");
+		console.info("Arguments: z=", z, " Boundaries: ", arguments[1], " NaN correction mode: ", nanCorrectionMode);
 		
-		let wl = this.getWaterline(z, 0);
+		let wl = this.getWaterline(z, nanCorrectionMode);
 		console.info("wl: ", wl); //DEBUG
 
 		let LOA = this.attributes.LOA;
@@ -1008,8 +1060,8 @@ Object.assign(Hull.prototype, {
 			sts[i] *= LOA;
 		}
 		
-		let hasMinX = (minX !== undefined);
-		let hasMaxX = (maxX !== undefined);
+		let hasMinX = (minX !== undefined) && minX!==sts[0];
+		let hasMaxX = (maxX !== undefined) && maxX!==sts[sts.length-1];
 		if (hasMinX || hasMaxX) {
 			let first=0;
 			let wlpre;
@@ -1135,12 +1187,12 @@ Object.assign(Hull.prototype, {
 				Cv: {x:0, y:0, z:0}
 			}) {
 			
-			let wlc = hull.waterlineCalculation(z);
+			let wlc = hull.waterlineCalculation(z,{},3);
 			let lev = {};
 			Object.assign(lev, wlc);
 			//Projected area calculation (approximate):
 			lev.prMinY = wlc.minY || 0;
-			lev.prMaxY = wlc.maxY || 0; //|| 0 right?
+			lev.prMaxY = wlc.maxY || 0;
 			lev.Ap = prev.Ap
 				+ trapezoidCalculation(prev.prMinY, prev.prMaxY, lev.prMinY, lev.prMaxY, prev.z, lev.z)["A"];
 			
@@ -1158,10 +1210,9 @@ Object.assign(Hull.prototype, {
 			//Find bilinear patches in the slice, and combine them.
 			//Many possibilities for getting the coordinate systems wrong.
 			let calculations = [];
-			//let wls = hull.halfBreadths.waterlines.map(wl=>wl*hull.attributes.Depth);
 			let sts = hull.halfBreadths.stations.map(st=>st*hull.attributes.LOA);
-			let wl = hull.getWaterline(z,0);
-			let prwl = hull.getWaterline(prev.z,0);
+			let wl = hull.getWaterline(z,3);
+			let prwl = hull.getWaterline(prev.z,3);
 			for (let j = 0; j < sts.length-1; j++) {
 				let port = 
 					bilinearPatchColumnCalculation(sts[j], sts[j+1], prev.z, z, -prwl[j], -wl[j], -prwl[j+1], -wl[j+1]);
@@ -1191,9 +1242,11 @@ Object.assign(Hull.prototype, {
 		
 		//Here is the returned function calculateAttributesAtDraft(T):
 		return function(T) {
-			if (T === undefined) {
+			if (isNaN(T)) {
 				console.error("Hull.prototype.calculateAttributesAtDraft(T): No draft specified. Returning undefined.");
 				return;
+			} else if (T<0 || T>this.attributes.Depth) {
+				console.error("Hull.prototype.calculateAttributesAtDraft(T): Draft parameter " + T + "outside valid range of [0,Depth]. Returning undefined.");
 			}
 			
 			let wls = this.halfBreadths.waterlines.map(wl=>this.attributes.Depth*wl);
diff --git a/examples/blockCase_compare_json.html b/examples/blockCase_compare_json.html
index a21cf53..96dd71a 100644
--- a/examples/blockCase_compare_json.html
+++ b/examples/blockCase_compare_json.html
@@ -99,7 +99,6 @@
 
 	
 	//Simplified calculations for triangular prism case, without trim:
-	//Assuming that T/*=*/wls.length-1) {
+		} else if (a >= wls.length-1) {
 			if (nanCorrectionMode===0) {
 				console.warn("getWaterLine: z above highest defined waterline. Defaulting to zeros.");
 				return new Array(sts.length).fill(0);
@@ -86,7 +88,7 @@ Object.assign(Hull.prototype, {
 				console.warn("getWaterLine: z above highest defined waterline. Outputting NaNs.");
 				return new Array(sts.length).fill(null);
 			}
-			else /*nanCorrectionMode===2*/ {
+			else /*nanCorrectionMode===2 || nanCorrectionMode===3*/ {
 				console.warn("getWaterLine: z above highest defined waterline. Proceeding with highest data entry.");
 				a = wls.length-2; //if this level is defined...
 				mu=1;
@@ -99,7 +101,7 @@ Object.assign(Hull.prototype, {
 		for (let j = 0; j < wl.length; j++) {
 			if (nanCorrectionMode === 0) {
 				if (a+1 > wls.length-1) {
-					wl[j] = lerp(tab[a][j], 0, 0.5);
+					wl[j] = lerp(tab[a][j], 0, 0.5); //Suspicious!
 				} else {
 					wl[j] = lerp(tab[a][j] || 0, tab[a+1][j] || 0, mu || 0.5);
 				}
@@ -109,7 +111,7 @@ Object.assign(Hull.prototype, {
 				} else {
 					wl[j] = lerp(tab[a][j], tab[a+1][j], mu);
 				}
-			} else {
+			} else if (nanCorrectionMode === 2) {
 				//If necessary, sample from below
 				let b = a;
 				while (b>0 && isNaN(tab[b][j])) {
@@ -141,14 +143,62 @@ Object.assign(Hull.prototype, {
 						upper = tab[c][j];
 					}
 				}
-				mu = c===b ? 0 : (a+(mu||0.5)-b)/(c-b);
+				mu = c===b ? 0 : (a+(mu||0.5)-b)/(c-b); //what is this? Some kind of self-explanatory code, perhaps?
 				wl[j] = lerp(lower, upper, mu);
+			} else /*nanCorrectionMode === 3*/ {
+				let lower, upper;
+				let b;
+				if (!isNaN(tab[a][j])) {
+					lower = tab[a][j];					
+				} else {
+					b = a+1;
+					while(b < wls.length && isNaN(tab[b][j])) {
+						b++;
+					}
+					if (b !== wls.length) {
+						//Inner NaN
+						lower = 0;
+					} else {
+						//Upper NaN, search below:
+						b = a-1;
+						while (b >= 0 && isNaN(tab[b][j])) {
+							b--;
+						}
+						if (b===-1) {
+							//No number found:
+							lower = 0;
+							upper = 0;
+						} else {
+							lower = tab[b][j];
+							upper = lower;
+						}
+					}
+				}
+				if (upper !== undefined) {}
+				else if (!isNaN(tab[a+1][j])) {
+					upper = tab[a+1][j];
+				} else {
+					//The cell value is NaN.
+					//Upper is not defined.
+					//That means either tab[a][j] is a number
+					//or tab[a][j] is an inner NaN and
+					//there exists at least one number above it.
+					//In both cases I have to check above a+1.
+					b = a+2;
+					while (b < wls.length && isNaN(tab[b][j])) {
+						b++;
+					}
+					if (b === wls.length) upper = lower;
+					else {
+						upper = tab[b][j];
+					}
+				}
+				wl[j] = lerp(lower, upper, mu || 0.5);
 			}
 			
 			//Scale numerical values
 			if (!isNaN(wl[j])) wl[j] *= 0.5*ha.BOA;
 		}
-
 		return wl;
 	},
 	getStation: function(x) {
@@ -184,14 +234,14 @@ Object.assign(Hull.prototype, {
 	
 	//THIS is a candidate for causing wrong Ix, Iy values.
 	//Much logic that can go wrong.
-										//typically deck bounds
-	waterlineCalculation: function(z, bounds) {
+									//typically deck bounds
+	waterlineCalculation: function(z, bounds, nanCorrectionMode=3) {
 		let {minX, maxX, minY, maxY} = bounds || {};
 
-		console.groupCollapsed("waterlineCalculation.");
-		console.info("Arguments: z=", z, " Boundaries: ", arguments[1]);
+		console.group/*Collapsed*/("waterlineCalculation.");
+		console.info("Arguments: z=", z, " Boundaries: ", arguments[1], " NaN correction mode: ", nanCorrectionMode);
 		
-		let wl = this.getWaterline(z, 0);
+		let wl = this.getWaterline(z, nanCorrectionMode);
 		console.info("wl: ", wl); //DEBUG
 
 		let LOA = this.attributes.LOA;
@@ -201,8 +251,8 @@ Object.assign(Hull.prototype, {
 			sts[i] *= LOA;
 		}
 		
-		let hasMinX = (minX !== undefined);
-		let hasMaxX = (maxX !== undefined);
+		let hasMinX = (minX !== undefined) && minX!==sts[0];
+		let hasMaxX = (maxX !== undefined) && maxX!==sts[sts.length-1];
 		if (hasMinX || hasMaxX) {
 			let first=0;
 			let wlpre;
@@ -328,12 +378,12 @@ Object.assign(Hull.prototype, {
 				Cv: {x:0, y:0, z:0}
 			}) {
 			
-			let wlc = hull.waterlineCalculation(z);
+			let wlc = hull.waterlineCalculation(z,{},3);
 			let lev = {};
 			Object.assign(lev, wlc);
 			//Projected area calculation (approximate):
 			lev.prMinY = wlc.minY || 0;
-			lev.prMaxY = wlc.maxY || 0; //|| 0 right?
+			lev.prMaxY = wlc.maxY || 0;
 			lev.Ap = prev.Ap
 				+ trapezoidCalculation(prev.prMinY, prev.prMaxY, lev.prMinY, lev.prMaxY, prev.z, lev.z)["A"];
 			
@@ -351,10 +401,9 @@ Object.assign(Hull.prototype, {
 			//Find bilinear patches in the slice, and combine them.
 			//Many possibilities for getting the coordinate systems wrong.
 			let calculations = [];
-			//let wls = hull.halfBreadths.waterlines.map(wl=>wl*hull.attributes.Depth);
 			let sts = hull.halfBreadths.stations.map(st=>st*hull.attributes.LOA);
-			let wl = hull.getWaterline(z,0);
-			let prwl = hull.getWaterline(prev.z,0);
+			let wl = hull.getWaterline(z,3);
+			let prwl = hull.getWaterline(prev.z,3);
 			for (let j = 0; j < sts.length-1; j++) {
 				let port = 
 					bilinearPatchColumnCalculation(sts[j], sts[j+1], prev.z, z, -prwl[j], -wl[j], -prwl[j+1], -wl[j+1]);
@@ -384,9 +433,11 @@ Object.assign(Hull.prototype, {
 		
 		//Here is the returned function calculateAttributesAtDraft(T):
 		return function(T) {
-			if (T === undefined) {
+			if (isNaN(T)) {
 				console.error("Hull.prototype.calculateAttributesAtDraft(T): No draft specified. Returning undefined.");
 				return;
+			} else if (T<0 || T>this.attributes.Depth) {
+				console.error("Hull.prototype.calculateAttributesAtDraft(T): Draft parameter " + T + "outside valid range of [0,Depth]. Returning undefined.");
 			}
 			
 			let wls = this.halfBreadths.waterlines.map(wl=>this.attributes.Depth*wl);
diff --git a/source/CoreClasses/Ship.js b/source/CoreClasses/Ship.js
index 0414559..b8f1131 100644
--- a/source/CoreClasses/Ship.js
+++ b/source/CoreClasses/Ship.js
@@ -17,7 +17,7 @@ Object.assign(Ship.prototype, {
 	constructor: Ship,
 	setFromSpecification: function(specification) {
 		this.attributes = specification.attributes || {};
-		this.structure = new Structure(specification.structure,this);
+		this.structure = new Structure(specification.structure/*,this*/);
 		//baseObjects and derivedObjects are arrays in the specification, but are objects (hashmaps) in the constructed ship object:
 		this.baseObjects = {};
 		for (let i = 0; i < specification.baseObjects.length; i++) {
diff --git a/source/CoreClasses/Structure.js b/source/CoreClasses/Structure.js
index 3f42d2d..26b677b 100644
--- a/source/CoreClasses/Structure.js
+++ b/source/CoreClasses/Structure.js
@@ -1,7 +1,7 @@
 //@EliasHasle
 
-function Structure(spec, ship) {
-	this.ship = ship;
+function Structure(spec/*, ship*/) {
+	//this.ship = ship;
 	JSONSpecObject.call(this, spec);
 }
 Structure.prototype = Object.create(JSONSpecObject.prototype);
@@ -63,7 +63,7 @@ Object.assign(Structure.prototype, {
 			let zc = d.zFloor+0.5*d.thickness;
 			let yc = d.yCentre;
 			let b = d.breadth;
-			let wlc = this.hull.waterlineCalculation(zc, {minX: d.xAft, maxX: d.xFwd, minY: yc-0.5*b, maxY: yc+0.5*b});
+			let wlc = this.hull.waterlineCalculation(zc, {minX: d.xAft, maxX: d.xFwd, minY: yc-0.5*b, maxY: yc+0.5*b}, 3);
 			components.push({
 				//approximation
 				mass: wlc.Awp*d.thickness*d.density,
diff --git a/source/functions/areaCalculations.js b/source/functions/areaCalculations.js
index 5f7db2f..24f4373 100644
--- a/source/functions/areaCalculations.js
+++ b/source/functions/areaCalculations.js
@@ -2,7 +2,7 @@
 
 //All inputs are numbers. The axes are given by a single coordinate.
 function steiner(I, A, sourceAxis, targetAxis) {
-	return I + A*(sourceAxis-targetAxis)^2;
+	return I + A*(sourceAxis-targetAxis)**2;
 }
 
 //Calculate area, center, Ix, Iy.
@@ -17,18 +17,18 @@ function trapezoidCalculation(xbase0, xbase1, xtop0, xtop1, ybase, ytop) {
 	let yc = (a==0 && b==0) ? ybase+0.5*h : ybase + h*(2*a+b)/(3*(a+b));
 	let d = xbase0+0.5*a; //shorthand
 	let xc = h===0 ? 0.25*(xbase0+xbase1+xtop0+xtop1) : d + (xtop0+0.5*b-d)*(yc-ybase)/h;
-	let Ix = (a==0 && b== 0) ? 0 : h^3*(a^2+4*a*b+b^2)/(36*(a+b));
+	let Ix = (a==0 && b== 0) ? 0 : h**3*(a**2+4*a*b+b**2)/(36*(a+b));
 
 	//For Iy I must decompose (I think negative results will work fine):
 	let Art1 = 0.5*(xtop0-xbase0)*h;
 	let xcrt1 = xbase0 + (xtop0-xbase0)/3;
-	let Iyrt1 = (xtop0-xbase0)^3*h/36;
+	let Iyrt1 = (xtop0-xbase0)**3*h/36;
 	let Arec = (xbase1-xtop0)*h;
 	let xcrec = 0.5*(xtop0+xbase1);
-	let Iyrec = (xbase1-xtop0)^3*h/12;
+	let Iyrec = (xbase1-xtop0)**3*h/12;
 	let Art2 = 0.5*(xbase1-xtop1)*h;
 	let xcrt2 = (xtop1 + (xbase1-xtop1)/3);
-	let Iyrt2 = (xbase1-xtop1)^3*h/36;
+	let Iyrt2 = (xbase1-xtop1)**3*h/36;
 
 	let Iy = steiner(Iyrt1, Art1, xcrt1, xc)
 		+ steiner(Iyrec, Arec, xcrec, xc)
@@ -85,10 +85,8 @@ function combineAreas(array) {
 
 //x and y here refers to coordinates in the plane that is being calculated on.
 function sectionCalculation({xs, ymins, ymaxs}) {
-	console.groupCollapsed("sectionCalculation");
+	console.group/*Collapsed*/("sectionCalculation");
 	console.info("Arguments (xs, ymins, ymaxs): ", arguments[0]);
-	
-	//Needed for Cwp (not a very efficient calculation, maybe):
 
 	let calculations = [];
 	for (let i = 0; i < xs.length-1; i++) {
diff --git a/source/functions/interpolation.js b/source/functions/interpolation.js
index ed95a78..3101738 100644
--- a/source/functions/interpolation.js
+++ b/source/functions/interpolation.js
@@ -2,22 +2,6 @@
 
 //Some interpolation helpers. Only linear and bilinear for now.
 
-//linear interpolation
-//Defaults are not finally decided
-//returns NaN if a and b are NaN or mu is NaN.
-function lerp(a, b, mu=0.5) {
-	if (isNaN(a)) return b;
-	if (isNaN(b)) return a;
-	return (1-mu)*a+mu*b;
-}
-
-//Test. Not safe yet.
-function linearFromArrays(xx, yy, x) {
-	let {index, mu} = bisectionSearch(xx, x);
-	if (index === undefined || mu === undefined) return 0;
-	return lerp(yy[index], yy[index+1], mu);	
-}
-
 /*Function that takes a sorted array as input, and finds the last index that holds a numerical value less than, or equal to, a given value.
 Returns an object with the index and an interpolation parameter mu that gives the position of value between index and index+1.
 */
@@ -44,6 +28,22 @@ function bisectionSearch(array, value) {
 	return {index, mu};
 }
 
+//linear interpolation
+//Defaults are not finally decided
+//returns NaN if a and b are NaN or mu is NaN.
+function lerp(a, b, mu=0.5) {
+	if (isNaN(a)) return b;
+	if (isNaN(b)) return a;
+	return (1-mu)*a+mu*b;
+}
+
+//Test. Not safe yet.
+function linearFromArrays(xx, yy, x) {
+	let {index, mu} = bisectionSearch(xx, x);
+	if (index === undefined || mu === undefined) return 0;
+	return lerp(yy[index], yy[index+1], mu);	
+}
+
 function bilinearUnitSquareCoeffs(z00, z01, z10, z11) {
 	let a00 = z00;
 	let a10 = z10-z00;
diff --git a/source/functions/volumeCalculations.js b/source/functions/volumeCalculations.js
index 30c906b..837a0e7 100644
--- a/source/functions/volumeCalculations.js
+++ b/source/functions/volumeCalculations.js
@@ -11,19 +11,19 @@ function bilinearPatchColumnCalculation(x1, x2, y1, y2, z11, z12, z21, z22) {
 	= INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x) dx]
 	= a00 + 0.5*a10 + 0.5*a01 + 0.25*a11
 	*/
-	let Ab = X*Y;
+	let Ab = X*Y; //area of base of patch column
 	let zAvg = (a00 + 0.5*a10 + 0.5*a01 + 0.25*a11);
 	let V = Math.abs(Ab*zAvg); //new: absolute value
 	let zc = 0.5*zAvg;
 	/*
-	To find xc, I need to integrate x*z over the unit square, and scale and translate to world coordinates afterwards:
+	To find xc, I need to integrate x*z over the unit square, and scale and translate to ship coordinates afterwards:
 	INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x)*x dx]
 	= 0.5*a00 + a10/3 + 0.25*a01 + a11/6
 	Scale and translate:*/
-	let xc = y1 + X*(0.5*a00 + a10/3 + 0.25*a01 + a11/6)
+	let xc = x1 + X*(0.5*a00 + a10/3 + 0.25*a01 + a11/6);
 	
 	//Similar for yc:
-	let yc = y1 + Y*(0.5*a00 + 0.25*a10 + a01/3 + a11/6)
+	let yc = y1 + Y*(0.5*a00 + 0.25*a10 + a01/3 + a11/6);
 	
 	//new: absolute value (OK?)
 	let As = Math.abs(bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22));
@@ -63,6 +63,8 @@ function combineVolumes(array) {
 }
 
 //For wetted area. I think this is right, but it is not tested.
+//The numerical integral will not scale well with larger geometries.
+//Then the full analytical solution is needed.
 function bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22, segs=10) {
 	let [b00,b10,b01,b11] = bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22);
 	/*
@@ -76,6 +78,8 @@ function bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22, segs=10) {
 	Tx X Ty = (-(b10+b11*y), -(b01+b11*x), 1)
 	|Tx X Ty| = sqrt((b10+b11*y)^2 + (b01+b11*x)^2 + 1)
 	
+	Now, to get the area I need to integrate |Tx X Ty| over X,Y.
+	
 	Wolfram Alpha gave me this for the inner integral using x (indefinite):
 	integral sqrt((b01 + b11 x)^2 + 1 + (b10+b11*y)^2) dx = ((b01 + b11*x) sqrt((b01 + b11*x)^2 + 1 + (b10+b11*y)^2) + (1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x))/(2*b11) + constant
 	That means this for the definite integral:

From 13b748fb69399dc804967101c875a984db497402 Mon Sep 17 00:00:00 2001
From: EliasHasle 
Date: Fri, 27 Oct 2017 15:22:06 +0200
Subject: [PATCH 11/20] Many changes, not done, don't apply to master branch

---
 build/Vessel.js                        | 194 +++++++++++++++----------
 examples/index.html                    |   2 +-
 source/CoreClasses/Hull.js             |  29 ++--
 source/functions/areaCalculations.js   |  72 +++++++++
 source/functions/volumeCalculations.js | 160 ++++++++------------
 5 files changed, 270 insertions(+), 187 deletions(-)

diff --git a/build/Vessel.js b/build/Vessel.js
index a9efb52..e638748 100644
--- a/build/Vessel.js
+++ b/build/Vessel.js
@@ -1,4 +1,4 @@
-//Vessel.js library, built 2017-10-26 22:24:21.268970, Checksum: 69efcb4dadd22b34bfc962e870f528a7
+//Vessel.js library, built 2017-10-27 15:17:23.678874, Checksum: d15ff7f6a612a6565c1b1e1ff341031a
 /*
 Import like this in HTML:
 
@@ -263,68 +263,6 @@ function sectionCalculation({xs, ymins, ymaxs}) {
 	console.info("Output: ", output);
 	console.groupEnd();
 	return output;
-}//@EliasHasle
-
-function bilinearPatchColumnCalculation(x1, x2, y1, y2, z11, z12, z21, z22) {
-	let X = x2-x1;
-	let Y = y2-y1;
-	let [a00, a10, a01, a11] = bilinearUnitSquareCoeffs(z11, z12, z21, z22);
-	/*
-	From here I call mux for x, and muy for y.
-	Integral over unit square:
-	INT[x from 0 to 1, INT[y from 0 to 1, (a00 + a10*x + a01*y + a11*x*y) dy] dx]
-	= INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x) dx]
-	= a00 + 0.5*a10 + 0.5*a01 + 0.25*a11
-	*/
-	let Ab = X*Y; //area of base of patch column
-	let zAvg = (a00 + 0.5*a10 + 0.5*a01 + 0.25*a11);
-	let V = Math.abs(Ab*zAvg); //new: absolute value
-	let zc = 0.5*zAvg;
-	/*
-	To find xc, I need to integrate x*z over the unit square, and scale and translate to ship coordinates afterwards:
-	INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x)*x dx]
-	= 0.5*a00 + a10/3 + 0.25*a01 + a11/6
-	Scale and translate:*/
-	let xc = x1 + X*(0.5*a00 + a10/3 + 0.25*a01 + a11/6);
-	
-	//Similar for yc:
-	let yc = y1 + Y*(0.5*a00 + 0.25*a10 + a01/3 + a11/6);
-	
-	//new: absolute value (OK?)
-	let As = Math.abs(bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22));
-	
-	return {Ab: Ab, As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
-}
-
-//Input: array of objects with calculation results for elements.
-//Output: the combined results.
-function combineVolumes(array) {
-	let V = 0;
-	let As = 0;
-	let Cv = {x:0, y:0, z:0};
-	let L = array.length;
-	if (L===0) return {V,As,Cv};
-	for (let i = 0; i < L; i++) {
-		let e = array[i];
-		V += e.V;
-		As += e.As; //typically wetted area
-		Cv.x += e.Cv.x*e.V;
-		Cv.y += e.Cv.y*e.V;
-		Cv.z += e.Cv.z*e.V;
-	}
-	//Safe zero check?
-	if (V!==0) {
-		Cv.x /= V;
-		Cv.y /= V;
-		Cv.z /= V;
-	} else {
-		console.warn("Zero volume combination.");
-		Cv.x /= L;
-		Cv.y /= L;
-		Cv.z /= L;
-	}
-	
-	return {V,As,Cv};//{V: V, As: As, Cv: Cv};
 }
 
 //For wetted area. I think this is right, but it is not tested.
@@ -383,7 +321,106 @@ function bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22, segs=10) {
 	A *= X*Y/(N*M); //dx dy
 	
 	return A;
-}//@MrEranwe
+}
+
+//Calculates the (arithmetic) average of the area of the two possible triangulations of the quad element (using two triangles).
+function elementArea(v1,v2,v3,v4) {
+	let A = 0.5*(Math.abs(signedTriangleArea(v1,v2,v3)) + Math.abs(signedTriangleArea(v2,v3,v4)) + Math.abs(signedTriangleArea(v3,v4,v1)) + Math.abs(signedTriangleArea(v4,v1,v2)));
+	return A;
+}
+
+function signedTriangleArea(v1,v2,v3) {
+	let u = addVec(v2,scaleVec(v1,-1));
+	let v = addVec(v3,scaleVec(v1,-1));
+	let c = crossProduct(u,v);
+	let A = 0.5*vecNorm(c);
+	return A;
+}//@EliasHasle
+
+//This one is broken, apparently.
+//It returns negative xc and yc when the z values are negative, 
+//even though all x and y values are positive.
+//I am currently doing some tests here of an (over-)simplified calculation
+//The results so far indicate that, for the prism hull, the results are almost identical, except that with the simple calculation the center of volume is almost right (but wrong enough to disqualify such a simple calculation). The bug causing very wrong (too big) submerged volume is likely to be found elsewhere.
+											  // xy
+function patchColumnCalculation(x1, x2, y1, y2, z11, z12, z21, z22) {
+	//DEBUG START:
+	//Simpler approximate calculation of volume:
+	let z = 0.25*(z11+z12+z21+z22);
+	let V = Math.abs((x2-x1)*(y2-y1)*z);
+
+	//Very approximate center of volume
+	//(does not account for different signs on z values,
+	//but that should be OK for hull offsets)
+	let xc = (x1*(z11+z12)+x2*(z21+z22))/(z11+z12+z21+z22);
+	let yc = (y1*(z11+z21)+y2*(z12+z22))/(z11+z12+z21+z22);
+	let zc = 0.5*z;
+	
+	//Simple triangle average approximation for area
+	let As = elementArea(
+		{x: x1, y: y1, z: z11},
+		{x: x1, y: y2, z: z12},
+		{x: x2, y: y1, z: z21},
+		{x: x2, y: y2, z: z22});
+	
+	return {As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
+	
+	//DEBUG END
+	
+	//Calculation based on a bilinear patch:
+
+	// let X = x2-x1;
+	// let Y = y2-y1;
+	// let [a00, a10, a01, a11] = bilinearUnitSquareCoeffs(z11, z12, z21, z22);
+	// /*
+	// From here I call mux for x, and muy for y.
+	// Integral over unit square:
+	// INT[x from 0 to 1, INT[y from 0 to 1, (a00 + a10*x + a01*y + a11*x*y) dy] dx]
+	// = INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x) dx]
+	// = a00 + 0.5*a10 + 0.5*a01 + 0.25*a11
+	// */
+	// let Ab = X*Y; //area of base of patch column
+	// let zAvg = (a00 + 0.5*a10 + 0.5*a01 + 0.25*a11);
+	// let V = Math.abs(Ab*zAvg); //new: absolute value
+	// let zc = 0.5*zAvg;
+	// /*
+	// To find xc, I need to integrate x*z over the unit square, and scale and translate to ship coordinates afterwards:
+	// INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x)*x dx]
+	// = 0.5*a00 + a10/3 + 0.25*a01 + a11/6
+	// Scale and translate:*/
+	// let xc = x1 + X*(0.5*a00 + a10/3 + 0.25*a01 + a11/6);
+	
+	// //Similar for yc:
+	// let yc = y1 + Y*(0.5*a00 + 0.25*a10 + a01/3 + a11/6);
+	
+	//new: absolute value (OK?)
+	//let As = Math.abs(bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22));
+	
+	// return {Ab: Ab, As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
+}
+
+//Input: array of objects with calculation results for elements.
+//Output: the combined results.
+function combineVolumes(array) {
+	let V = 0;
+	let As = 0;
+	let Cv = {x:0, y:0, z:0};
+	let L = array.length;
+	//if (L===0) return {V,As,Cv};
+	for (let i = 0; i < L; i++) {
+		let e = array[i];
+		V += e.V;
+		As += e.As; //typically wetted area
+		//console.log(e.Cv);
+		Cv = addVec(Cv, scaleVec(e.Cv, e.V));
+	}
+	Cv = scaleVec(Cv, 1/(V || L || 1));
+	
+	//console.info("combineVolumes: Combined Cv is (" + Cv.x + ", " + Cv.y + ", " + Cv.z + ").");
+	
+	return {V,As,Cv};//{V: V, As: As, Cv: Cv};
+}
+//@MrEranwe
 //@EliasHasle
 
 "use strict";
@@ -1215,25 +1252,28 @@ Object.assign(Hull.prototype, {
 			let prwl = hull.getWaterline(prev.z,3);
 			for (let j = 0; j < sts.length-1; j++) {
 				let port = 
-					bilinearPatchColumnCalculation(sts[j], sts[j+1], prev.z, z, -prwl[j], -wl[j], -prwl[j+1], -wl[j+1]);
+					patchColumnCalculation(sts[j], sts[j+1], prev.z, z, -prwl[j], -wl[j], -prwl[j+1], -wl[j+1]);
 				calculations.push(port);
 				let star =
-					bilinearPatchColumnCalculation(sts[j], sts[j+1], prev.z, z, prwl[j], wl[j], prwl[j+1], wl[j+1]);
+					patchColumnCalculation(sts[j], sts[j+1], prev.z, z, prwl[j], wl[j], prwl[j+1], wl[j+1]);
 				calculations.push(star);
 			}
 			let C = combineVolumes(calculations);
+			//Cv of slice. Note that switching of yz must
+			//be done before combining with previous level
+			let Cv = {x: C.Cv.x, y: C.Cv.z, z: C.Cv.y};
+			
 			lev.Vs = prev.Vs + C.V; //hull volume below z
 			lev.As = prev.As + C.As; //outside surface below z
 
-			//center of volume below z (some potential for accumulated rounding error):
-			let Cv = addVec(scaleVec(prev.Cv,prev.Vs),
-					scaleVec(C.Cv,C.V));
-			let V = prev.Vs+C.V;
-			if (V!==0) {
-				Cv = scaleVec(Cv, 1/(prev.Vs+C.V));
-			}
-						//Note switching of yz
-			lev.Cv = {x: Cv.x, y: Cv.z, z: Cv.y};
+			//Simple approximation (wrong) for end caps:
+			lev.As += 2*lev.Ap;
+			
+			//center of volume below z (some potential for accumulated rounding error when calculating an accumulated average like this):
+			lev.Cv = scaleVec(addVec(
+						scaleVec(prev.Cv,prev.Vs),
+						scaleVec(Cv,C.V)
+					), 1/(lev.Vs || 2));
 			
 			lev.Cb = lev.Vs/lev.Vbb;
 			
@@ -1265,7 +1305,9 @@ Object.assign(Hull.prototype, {
 			//Find highest data waterline below water:
 			let {index: previ} = bisectionSearch(wls, T);
 			
-			let lc = levelCalculation(this, T, this.levels[previ]);
+			//console.info("Highest data waterline below water: " + previ);
+			
+			let lc = levelCalculation(this, T, this.levels[previ] || undefined);
 			
 			//Filter and rename for output
 			return {
diff --git a/examples/index.html b/examples/index.html
index e977d40..9a51796 100644
--- a/examples/index.html
+++ b/examples/index.html
@@ -30,7 +30,7 @@ 
    Examples
    
 
    Ship3D_load_state_controls.html
     
    Very simple test of changing tank fullness with a slider to see change of draft.
    
 
    Vessel.js_test.html
     
    A script that does a few calculations based on the provided specification, and outputs to the screen and the console.
    
 
    blockCase_compare.html
     
    Comparison between manual and automatic calculations on a very simple ship.
    
-
    blockCase_compare_json.html
     
    Sketched alternative implementation of same comparison.
    
+
    blockCase_compare_json.html
     
    Sketched alternative implementation of same comparison.
    
 
    PX121_calculations.html
     
    Some calculations on a test specification partially based on an Ulstein vessel.
    
 
 
diff --git a/source/CoreClasses/Hull.js b/source/CoreClasses/Hull.js
index 93e4d64..a6ab1a5 100644
--- a/source/CoreClasses/Hull.js
+++ b/source/CoreClasses/Hull.js
@@ -406,25 +406,28 @@ Object.assign(Hull.prototype, {
 			let prwl = hull.getWaterline(prev.z,3);
 			for (let j = 0; j < sts.length-1; j++) {
 				let port = 
-					bilinearPatchColumnCalculation(sts[j], sts[j+1], prev.z, z, -prwl[j], -wl[j], -prwl[j+1], -wl[j+1]);
+					patchColumnCalculation(sts[j], sts[j+1], prev.z, z, -prwl[j], -wl[j], -prwl[j+1], -wl[j+1]);
 				calculations.push(port);
 				let star =
-					bilinearPatchColumnCalculation(sts[j], sts[j+1], prev.z, z, prwl[j], wl[j], prwl[j+1], wl[j+1]);
+					patchColumnCalculation(sts[j], sts[j+1], prev.z, z, prwl[j], wl[j], prwl[j+1], wl[j+1]);
 				calculations.push(star);
 			}
 			let C = combineVolumes(calculations);
+			//Cv of slice. Note that switching of yz must
+			//be done before combining with previous level
+			let Cv = {x: C.Cv.x, y: C.Cv.z, z: C.Cv.y};
+			
 			lev.Vs = prev.Vs + C.V; //hull volume below z
 			lev.As = prev.As + C.As; //outside surface below z
 
-			//center of volume below z (some potential for accumulated rounding error):
-			let Cv = addVec(scaleVec(prev.Cv,prev.Vs),
-					scaleVec(C.Cv,C.V));
-			let V = prev.Vs+C.V;
-			if (V!==0) {
-				Cv = scaleVec(Cv, 1/(prev.Vs+C.V));
-			}
-						//Note switching of yz
-			lev.Cv = {x: Cv.x, y: Cv.z, z: Cv.y};
+			//Simple (and wrong) approximation for end caps:
+			lev.As += 2*lev.Ap; //(will work for prism case)
+			
+			//center of volume below z (some potential for accumulated rounding error when calculating an accumulated average like this):
+			lev.Cv = scaleVec(addVec(
+						scaleVec(prev.Cv,prev.Vs),
+						scaleVec(Cv,C.V)
+					), 1/(lev.Vs || 2));
 			
 			lev.Cb = lev.Vs/lev.Vbb;
 			
@@ -456,7 +459,9 @@ Object.assign(Hull.prototype, {
 			//Find highest data waterline below water:
 			let {index: previ} = bisectionSearch(wls, T);
 			
-			let lc = levelCalculation(this, T, this.levels[previ]);
+			//console.info("Highest data waterline below water: " + previ);
+			
+			let lc = levelCalculation(this, T, this.levels[previ] || undefined);
 			
 			//Filter and rename for output
 			return {
diff --git a/source/functions/areaCalculations.js b/source/functions/areaCalculations.js
index 24f4373..cda838e 100644
--- a/source/functions/areaCalculations.js
+++ b/source/functions/areaCalculations.js
@@ -106,4 +106,76 @@ function sectionCalculation({xs, ymins, ymaxs}) {
 	console.info("Output: ", output);
 	console.groupEnd();
 	return output;
+}
+
+//For wetted area. I think this is right, but it is not tested.
+//The numerical integral will not scale well with larger geometries.
+//Then the full analytical solution is needed.
+function bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22, segs=10) {
+	let [b00,b10,b01,b11] = bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22);
+	/*
+	z(x,y) = b00 + b10*x + b01*y + b11*xy
+	Partial derivative in x: b10 + b11*y
+	Partial derivative in y: b01 + b11*x
+	I think this will be right:
+	Tx(y) = (1, 0, b10+b11*y)
+	Ty(x) = (0, 1, b01+b11*x)
+	Then:
+	Tx X Ty = (-(b10+b11*y), -(b01+b11*x), 1)
+	|Tx X Ty| = sqrt((b10+b11*y)^2 + (b01+b11*x)^2 + 1)
+	
+	Now, to get the area I need to integrate |Tx X Ty| over X,Y.
+	
+	Wolfram Alpha gave me this for the inner integral using x (indefinite):
+	integral sqrt((b01 + b11 x)^2 + 1 + (b10+b11*y)^2) dx = ((b01 + b11*x) sqrt((b01 + b11*x)^2 + 1 + (b10+b11*y)^2) + (1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x))/(2*b11) + constant
+	That means this for the definite integral:
+	((b01 + b11*x2)*sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2) + 1 + (b10+b11*y)^2*ln(sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x2))/(2*b11) - ((b01 + b11*x1)*sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2) + (1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x1))/(2*b11)
+	=
+	(b01 + b11*x2)*sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2)/(2*b11)
+	+(1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x2)/(2*b11))
+	- (b01 + b11*x1)*sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2)/(2*b11)
+	- (1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x1)/(2*b11)
+	=
+	(b01 + b11*x2)*sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2)/(2*b11)
+	- (b01 + b11*x1)*sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2)/(2*b11)
+	+(1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x2)/(2*b11))
+	- (1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x1)/(2*b11)
+	
+	The two first integrals are similar, and the two last are also similar. With A=+-(b01 + b11*xi)/(2*b11), B=(b01 + b11*xi)^2+1, C=b10 and D=b11 (where xi represents either x1 or x2, and +- represents + for x2 and - for x1), I can calculate the integral of sqrt(B+(C+D*y)^2) and multiply by A. That integral is on the same form as the first one.
+	
+	The two last integrals can be represented by setting A=+-1/(2*b11), B=(b01 + b11*xi)^2+1, C=b01+b11*xi, D=b10, E=b11, and calculating the integral of (1+(D+E*y)^2)*ln(sqrt(B+(D+E*y)^2)+C), and multiplying by A.
+	Here is the integration result from Wolfram Alpha:
+	integral(1 + (D + E y)^2) log(sqrt(B + (D + E y)^2) + C) dy = (-(6 (B^2 - 2 B C^2 - 3 B + C^4 + 3 C^2) tan^(-1)((D + E y)/sqrt(B - C^2)))/sqrt(B - C^2) + (6 (B^2 - 2 B C^2 - 3 B + C^4 + 3 C^2) tan^(-1)((C (D + E y))/(sqrt(B - C^2) sqrt(B + (D + E y)^2))))/sqrt(B - C^2) + 6 (B - C^2 - 3) (D + E y) + 3 C (-3 B + 2 C^2 + 6) log(sqrt(B + (D + E y)^2) + D + E y) + 3 C (D + E y) sqrt(B + (D + E y)^2) + 6 ((D + E y)^2 + 3) (D + E y) log(sqrt(B + (D + E y)^2) + C) - 2 (D + E y)^3)/(18 E) + constant
+
+	I am glad I did not try to do this by hand. But combining these formulae, we can get an exact integral of the area of a bilinear patch. Later. Bilinear patches are not an exact representation anyway. We may opt for something else.
+	*/
+	
+	//Simple numerical calculation of double integral:
+	let A = 0;
+	let X = x2-x1, Y = y2-y1;
+	let N = segs, M = segs;
+	for (let i = 0; i < N; i++) {
+		let x = x1 + ((i+0.5)/N)*X;
+		for (let j = 0; j < M; j++) {
+			let y = y1 + ((j+0.5)/M)*Y;
+			A += Math.sqrt((b10+b11*y)**2 + (b01+b11*x)**2 + 1);
+		}
+	}
+	A *= X*Y/(N*M); //dx dy
+	
+	return A;
+}
+
+//Calculates the (arithmetic) average of the area of the two possible triangulations of the quad element (using two triangles).
+function elementArea(v1,v2,v3,v4) {
+	let A = 0.5*(Math.abs(signedTriangleArea(v1,v2,v3)) + Math.abs(signedTriangleArea(v2,v3,v4)) + Math.abs(signedTriangleArea(v3,v4,v1)) + Math.abs(signedTriangleArea(v4,v1,v2)));
+	return A;
+}
+
+function signedTriangleArea(v1,v2,v3) {
+	let u = addVec(v2,scaleVec(v1,-1));
+	let v = addVec(v3,scaleVec(v1,-1));
+	let c = crossProduct(u,v);
+	let A = 0.5*vecNorm(c);
+	return A;
 }
\ No newline at end of file
diff --git a/source/functions/volumeCalculations.js b/source/functions/volumeCalculations.js
index 837a0e7..e036d33 100644
--- a/source/functions/volumeCalculations.js
+++ b/source/functions/volumeCalculations.js
@@ -1,34 +1,65 @@
 //@EliasHasle
 
-function bilinearPatchColumnCalculation(x1, x2, y1, y2, z11, z12, z21, z22) {
-	let X = x2-x1;
-	let Y = y2-y1;
-	let [a00, a10, a01, a11] = bilinearUnitSquareCoeffs(z11, z12, z21, z22);
-	/*
-	From here I call mux for x, and muy for y.
-	Integral over unit square:
-	INT[x from 0 to 1, INT[y from 0 to 1, (a00 + a10*x + a01*y + a11*x*y) dy] dx]
-	= INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x) dx]
-	= a00 + 0.5*a10 + 0.5*a01 + 0.25*a11
-	*/
-	let Ab = X*Y; //area of base of patch column
-	let zAvg = (a00 + 0.5*a10 + 0.5*a01 + 0.25*a11);
-	let V = Math.abs(Ab*zAvg); //new: absolute value
-	let zc = 0.5*zAvg;
-	/*
-	To find xc, I need to integrate x*z over the unit square, and scale and translate to ship coordinates afterwards:
-	INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x)*x dx]
-	= 0.5*a00 + a10/3 + 0.25*a01 + a11/6
-	Scale and translate:*/
-	let xc = x1 + X*(0.5*a00 + a10/3 + 0.25*a01 + a11/6);
+//This one is broken, apparently.
+//It returns negative xc and yc when the z values are negative, 
+//even though all x and y values are positive.
+//I am currently doing some tests here of an (over-)simplified calculation
+//The results so far indicate that, for the prism hull, the results are almost identical, except that with the simple calculation the center of volume is almost right (but wrong enough to disqualify such a simple calculation). The bug causing very wrong (too big) submerged volume is likely to be found elsewhere.
+											  // xy
+function patchColumnCalculation(x1, x2, y1, y2, z11, z12, z21, z22) {
+	//DEBUG START:
+	//Simpler approximate calculation of volume:
+	let z = 0.25*(z11+z12+z21+z22);
+	let V = Math.abs((x2-x1)*(y2-y1)*z);
+
+	//Very approximate center of volume
+	//(does not account for different signs on z values,
+	//but that should be OK for hull offsets)
+	let xc = (x1*(z11+z12)+x2*(z21+z22))/(z11+z12+z21+z22);
+	let yc = (y1*(z11+z21)+y2*(z12+z22))/(z11+z12+z21+z22);
+	let zc = 0.5*z;
+	
+	//Simple triangle average approximation for area
+	let As = elementArea(
+		{x: x1, y: y1, z: z11},
+		{x: x1, y: y2, z: z12},
+		{x: x2, y: y1, z: z21},
+		{x: x2, y: y2, z: z22});
+	
+	return {As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
 	
-	//Similar for yc:
-	let yc = y1 + Y*(0.5*a00 + 0.25*a10 + a01/3 + a11/6);
+	//DEBUG END
+	
+	//Calculation based on a bilinear patch:
+
+	// let X = x2-x1;
+	// let Y = y2-y1;
+	// let [a00, a10, a01, a11] = bilinearUnitSquareCoeffs(z11, z12, z21, z22);
+	// /*
+	// From here I call mux for x, and muy for y.
+	// Integral over unit square:
+	// INT[x from 0 to 1, INT[y from 0 to 1, (a00 + a10*x + a01*y + a11*x*y) dy] dx]
+	// = INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x) dx]
+	// = a00 + 0.5*a10 + 0.5*a01 + 0.25*a11
+	// */
+	// let Ab = X*Y; //area of base of patch column
+	// let zAvg = (a00 + 0.5*a10 + 0.5*a01 + 0.25*a11);
+	// let V = Math.abs(Ab*zAvg); //new: absolute value
+	// let zc = 0.5*zAvg;
+	// /*
+	// To find xc, I need to integrate x*z over the unit square, and scale and translate to ship coordinates afterwards:
+	// INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x)*x dx]
+	// = 0.5*a00 + a10/3 + 0.25*a01 + a11/6
+	// Scale and translate:*/
+	// let xc = x1 + X*(0.5*a00 + a10/3 + 0.25*a01 + a11/6);
+	
+	// //Similar for yc:
+	// let yc = y1 + Y*(0.5*a00 + 0.25*a10 + a01/3 + a11/6);
 	
 	//new: absolute value (OK?)
-	let As = Math.abs(bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22));
+	//let As = Math.abs(bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22));
 	
-	return {Ab: Ab, As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
+	// return {Ab: Ab, As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
 }
 
 //Input: array of objects with calculation results for elements.
@@ -38,84 +69,17 @@ function combineVolumes(array) {
 	let As = 0;
 	let Cv = {x:0, y:0, z:0};
 	let L = array.length;
-	if (L===0) return {V,As,Cv};
+	//if (L===0) return {V,As,Cv};
 	for (let i = 0; i < L; i++) {
 		let e = array[i];
 		V += e.V;
 		As += e.As; //typically wetted area
-		Cv.x += e.Cv.x*e.V;
-		Cv.y += e.Cv.y*e.V;
-		Cv.z += e.Cv.z*e.V;
-	}
-	//Safe zero check?
-	if (V!==0) {
-		Cv.x /= V;
-		Cv.y /= V;
-		Cv.z /= V;
-	} else {
-		console.warn("Zero volume combination.");
-		Cv.x /= L;
-		Cv.y /= L;
-		Cv.z /= L;
+		//console.log(e.Cv);
+		Cv = addVec(Cv, scaleVec(e.Cv, e.V));
 	}
+	Cv = scaleVec(Cv, 1/(V || L || 1));
+	
+	//console.info("combineVolumes: Combined Cv is (" + Cv.x + ", " + Cv.y + ", " + Cv.z + ").");
 	
 	return {V,As,Cv};//{V: V, As: As, Cv: Cv};
 }
-
-//For wetted area. I think this is right, but it is not tested.
-//The numerical integral will not scale well with larger geometries.
-//Then the full analytical solution is needed.
-function bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22, segs=10) {
-	let [b00,b10,b01,b11] = bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22);
-	/*
-	z(x,y) = b00 + b10*x + b01*y + b11*xy
-	Partial derivative in x: b10 + b11*y
-	Partial derivative in y: b01 + b11*x
-	I think this will be right:
-	Tx(y) = (1, 0, b10+b11*y)
-	Ty(x) = (0, 1, b01+b11*x)
-	Then:
-	Tx X Ty = (-(b10+b11*y), -(b01+b11*x), 1)
-	|Tx X Ty| = sqrt((b10+b11*y)^2 + (b01+b11*x)^2 + 1)
-	
-	Now, to get the area I need to integrate |Tx X Ty| over X,Y.
-	
-	Wolfram Alpha gave me this for the inner integral using x (indefinite):
-	integral sqrt((b01 + b11 x)^2 + 1 + (b10+b11*y)^2) dx = ((b01 + b11*x) sqrt((b01 + b11*x)^2 + 1 + (b10+b11*y)^2) + (1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x))/(2*b11) + constant
-	That means this for the definite integral:
-	((b01 + b11*x2)*sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2) + 1 + (b10+b11*y)^2*ln(sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x2))/(2*b11) - ((b01 + b11*x1)*sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2) + (1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x1))/(2*b11)
-	=
-	(b01 + b11*x2)*sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2)/(2*b11)
-	+(1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x2)/(2*b11))
-	- (b01 + b11*x1)*sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2)/(2*b11)
-	- (1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x1)/(2*b11)
-	=
-	(b01 + b11*x2)*sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2)/(2*b11)
-	- (b01 + b11*x1)*sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2)/(2*b11)
-	+(1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x2)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x2)/(2*b11))
-	- (1 + (b10+b11*y)^2)*ln(sqrt((b01 + b11*x1)^2 + 1 + (b10+b11*y)^2) + b01 + b11*x1)/(2*b11)
-	
-	The two first integrals are similar, and the two last are also similar. With A=+-(b01 + b11*xi)/(2*b11), B=(b01 + b11*xi)^2+1, C=b10 and D=b11 (where xi represents either x1 or x2, and +- represents + for x2 and - for x1), I can calculate the integral of sqrt(B+(C+D*y)^2) and multiply by A. That integral is on the same form as the first one.
-	
-	The two last integrals can be represented by setting A=+-1/(2*b11), B=(b01 + b11*xi)^2+1, C=b01+b11*xi, D=b10, E=b11, and calculating the integral of (1+(D+E*y)^2)*ln(sqrt(B+(D+E*y)^2)+C), and multiplying by A.
-	Here is the integration result from Wolfram Alpha:
-	integral(1 + (D + E y)^2) log(sqrt(B + (D + E y)^2) + C) dy = (-(6 (B^2 - 2 B C^2 - 3 B + C^4 + 3 C^2) tan^(-1)((D + E y)/sqrt(B - C^2)))/sqrt(B - C^2) + (6 (B^2 - 2 B C^2 - 3 B + C^4 + 3 C^2) tan^(-1)((C (D + E y))/(sqrt(B - C^2) sqrt(B + (D + E y)^2))))/sqrt(B - C^2) + 6 (B - C^2 - 3) (D + E y) + 3 C (-3 B + 2 C^2 + 6) log(sqrt(B + (D + E y)^2) + D + E y) + 3 C (D + E y) sqrt(B + (D + E y)^2) + 6 ((D + E y)^2 + 3) (D + E y) log(sqrt(B + (D + E y)^2) + C) - 2 (D + E y)^3)/(18 E) + constant
-
-	I am glad I did not try to do this by hand. But combining these formulae, we can get an exact integral of the area of a bilinear patch. Later. Bilinear patches are not an exact representation anyway. We may opt for something else.
-	*/
-	
-	//Simple numerical calculation of double integral:
-	let A = 0;
-	let X = x2-x1, Y = y2-y1;
-	let N = segs, M = segs;
-	for (let i = 0; i < N; i++) {
-		let x = x1 + ((i+0.5)/N)*X;
-		for (let j = 0; j < M; j++) {
-			let y = y1 + ((j+0.5)/M)*Y;
-			A += Math.sqrt((b10+b11*y)**2 + (b01+b11*x)**2 + 1);
-		}
-	}
-	A *= X*Y/(N*M); //dx dy
-	
-	return A;
-}
\ No newline at end of file

From e7fd72bcaf5663b6c408c0992593d2a4c5b1f248 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Vicente=20Alejandro=20Iv=C3=A1=C3=B1ez=20Encinas?=
 
Date: Fri, 3 Nov 2017 09:57:00 +0100
Subject: [PATCH 12/20] 1st change

---
 examples/blockCase_compare.html               | 31 +++++++++++++++----
 .../data/ship_specifications/blockCase.json   |  5 +--
 2 files changed, 28 insertions(+), 8 deletions(-)

diff --git a/examples/blockCase_compare.html b/examples/blockCase_compare.html
index 720e381..6ad8972 100644
--- a/examples/blockCase_compare.html
+++ b/examples/blockCase_compare.html
@@ -212,6 +212,18 @@ 
    Comparation
    
                            
                             
                         
+                         
+                           Lightweigth (kg) 
+                            
+                           
+                            
+                        
+                        
+                           Deadweight (kg) 
+                            
+                           
+                            
+                        
                         
                     
                 
@@ -344,11 +356,16 @@ 
    Comparation
    
             let KBc = a.KB.toFixed(3);
             let KGc = a.KG.toFixed(3);
             let BMc = a.BM.toFixed(3);
-            let draftc = ship.calculateDraft().toFixed(3);
-            let lwlc = ship.structure.hull.calculateAttributesAtDraft().LWL;
-            let bwlc = ship.structure.hull.calculateAttributesAtDraft().BWL;
-            let Awpc = ship.structure.hull.calculateAttributesAtDraft().Awp;
-            let cbc = ship.structure.hull.calculateAttributesAtDraft().Cb;
+            let draftc = ship.designState.calculationParameters.Draft_design;
+            
+            //baseObjects not sure why is not working
+            let lightweight = ship.structure.hull.attributes.lightweight + 15000;
+            let deadweight = 14000 + 56000;
+            
+            let lwlc = ship.structure.hull.calculateAttributesAtDraft(draftc).LWL;
+            let bwlc = ship.structure.hull.calculateAttributesAtDraft(draftc).BWL;
+            let Awpc = ship.structure.hull.calculateAttributesAtDraft(draftc).Awp;
+            let cbc = ship.structure.hull.calculateAttributesAtDraft(draftc).Cb;
             
             document.getElementById("KBc").innerHTML = KBc;
             document.getElementById("KGc").innerHTML = KGc;
@@ -358,6 +375,8 @@ 
    Comparation
    
             document.getElementById("bwlc").innerHTML = bwlc.toFixed(3);
             document.getElementById("Awpc").innerHTML = Awpc.toFixed(3);
             document.getElementById("cbc").innerHTML = cbc.toFixed(3);
+            document.getElementById("lightweight").innerHTML = lightweight;
+            document.getElementById("deadweight").innerHTML = deadweight;
             
             //% diff
             let KBd = (1-KB/KBc)*100;
@@ -379,7 +398,7 @@ 
    Comparation
    
             document.getElementById("awpd").innerHTML = awpd.toFixed(3);
             
             
-            console.log(ship.structure.hull.calculateAttributesAtDraft())
+            //console.log(ship.structure.hull.calculateAttributesAtDraft(draftc))
 		}
 		
 		function animate() {
diff --git a/examples/data/ship_specifications/blockCase.json b/examples/data/ship_specifications/blockCase.json
index ee5e741..3f201e1 100644
--- a/examples/data/ship_specifications/blockCase.json
+++ b/examples/data/ship_specifications/blockCase.json
@@ -26,7 +26,8 @@
 			"LOA": 20,
 			"BOA": 10,
 			"Depth": 4,
-			"APP": 0
+			"APP": 0,
+            "lightweight": 15000
 		},
 		"halfBreadths": {
 			"waterlines": [0, 1],
@@ -60,7 +61,7 @@
 	}
 },
 "baseObjects": 
-[{"id":"TankL2.1B6.13H4.3fdundefinedFtank1.stl","affiliations":{},"boxDimensions":{"length":4,"breadth":6,"height":5},"capabilities":{},"file3D":"tank1.stl","baseState":{"fullness":1},"weightInformation":{"contentDensity":1000,"volumeCapacity":120,"lightweight":553.539,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}},{"id":"TankL2.1B6.13H4.3fdundefinedFundefined","affiliations":{},"boxDimensions":{"length":5,"breadth":4,"height":2},"capabilities":{},"baseState":{"fullness":0.5},"weightInformation":{"contentDensity":1000,"volumeCapacity":20,"lightweight":553.539,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}},{"id":"TankL4.199999999999999B6.7H4.3fdundefinedFtank2.stl","affiliations":{},"boxDimensions":{"length":10,"breadth":4,"height":2},"capabilities":{},"file3D":"tank2.stl","baseState":{"fullness":0.5},"weightInformation":{"contentDensity":1000,"volumeCapacity":48,"lightweight":1210.0199999999998,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}}],
+[{"id":"TankL2.1B6.13H4.3fdundefinedFtank1.stl","affiliations":{},"boxDimensions":{"length":4,"breadth":6,"height":5},"capabilities":{},"file3D":"tank1.stl","baseState":{"fullness":0},"weightInformation":{"contentDensity":0,"volumeCapacity":120,"lightweight":15000,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}},{"id":"TankL2.1B6.13H4.3fdundefinedFundefined","affiliations":{},"boxDimensions":{"length":5,"breadth":4,"height":2},"capabilities":{},"baseState":{"fullness":0},"weightInformation":{"contentDensity":0,"volumeCapacity":20,"lightweight":14000,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}},{"id":"TankL4.199999999999999B6.7H4.3fdundefinedFtank2.stl","affiliations":{},"boxDimensions":{"length":10,"breadth":4,"height":2},"capabilities":{},"file3D":"tank2.stl","baseState":{"fullness":0},"weightInformation":{"contentDensity":0,"volumeCapacity":48,"lightweight":56000,"fullnessCGMapping":{"fullnesses":[0,0.25,0.5,0.75,1],"cgs":[[0,0,2.15],[0,0,0.6449999999999999],[0,0,1.2899999999999998],[0,0,1.72],[0,0,2.15]]}}}],
 
 "derivedObjects": [{"id":"Tank1","baseObject":"TankL2.1B6.13H4.3fdundefinedFtank1.stl","affiliations":{"Deck":"WheelHouseTop","SFI":"101"},"referenceState":{"xCentre":2,"yCentre":0,"zBase":4}},{"id":"Tank2","baseObject":"TankL2.1B6.13H4.3fdundefinedFundefined","affiliations":{"Deck":"BridgeDeck","SFI":"102"},"referenceState":{"xCentre":5,"yCentre":0,"zBase":2}},{"id":"Tank3","baseObject":"TankL4.199999999999999B6.7H4.3fdundefinedFtank2.stl","affiliations":{"Deck":"BridgeDeck","SFI":"103"},"referenceState":{"xCentre":14,"yCentre":0,"zBase":2}}]
 

From b11d06786e51da39f4686d19012e90adac81c320 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Vicente=20Alejandro=20Iv=C3=A1=C3=B1ez=20Encinas?=
 
Date: Fri, 3 Nov 2017 11:51:07 +0100
Subject: [PATCH 13/20] New values

Add:
-Draft as input
-Displacement as output
-Few more values from hand calculations

To do:

-getWeight from baseObjects
-cb from program
-add rest of values from program
-finish structural
---
 examples/blockCase_compare.html | 86 +++++++++++++++++++++++++--------
 1 file changed, 66 insertions(+), 20 deletions(-)

diff --git a/examples/blockCase_compare.html b/examples/blockCase_compare.html
index 6ad8972..a59685d 100644
--- a/examples/blockCase_compare.html
+++ b/examples/blockCase_compare.html
@@ -143,7 +143,7 @@ 
    Comparation
    
                         
                            Displacement (t) 
                             
-                            
+                             
                             
                         
                         
@@ -215,12 +215,42 @@ 
    Comparation
    
                          
                            Lightweigth (kg) 
                             
-                           
+                           
                             
                         
                         
                            Deadweight (kg) 
                             
+                           
+                            
+                        
+                        
+                           Trim (m) 
+                            
+                           
+                            
+                        
+                        
+                           Trim (º) 
+                            
+                           
+                            
+                        
+                        
+                           Heel (m) 
+                            
+                           
+                            
+                        
+                        
+                           Maximum moment (ton m) 
+                            
+                           
+                            
+                        
+                        
+                           Shear  (t) 
+                            
                            
                             
                         
@@ -319,22 +349,28 @@ 
    Comparation
    
             document.getElementById("GM2").innerHTML = GMc;
             
             //hand values
-            let disp = 102.4;
-            let draftms = 1.989;
-            let draftap = 1.98;
-            let draftfp = 1.999;
+            let disp = 100;
+            let draftms = 2.028;
+            let draftap = 2.082;
+            let draftfp = 1.974;
             let lwl = 20;
-            let bwl = 5.046;
-            let Awet = 137.989;
-            let Awp = 99.95;
+            let bwl = 5.206;
+            let Awet = 136.989;
+            let Awp = 98.95;
             let Cb = disp / (draftms * lwl * bwl * 1.025);
-            let LCBx = 10.065;
-            let LCFx = 10.032;
-            let KB = 1.333;
-            let KG = 3.513;
-            let BM = 2.082;
+            let LCBx = 9.636;
+            let LCFx = 9.818;
+            let KB = 1.319;
+            let KG = 3.525;
+            let BM = 2.061;
             let GM = KB + BM - KG;
-            let trim = -0.019;
+            let trimd = 0.108;
+            let trima = 0.618;
+            let heel = 0;
+            let maxmom = 40.072;
+            let shear = 13.327;
+            let lightweight = 30000;
+            let deadweight = 70000;
             
             document.getElementById("disp").innerHTML = disp.toFixed(3);
             document.getElementById("draftms").innerHTML = draftms.toFixed(3);
@@ -343,7 +379,7 @@ 
    Comparation
    
             document.getElementById("lwl").innerHTML = lwl.toFixed(3);
             document.getElementById("bwl").innerHTML = bwl.toFixed(3);
             document.getElementById("Awet").innerHTML = Awet.toFixed(3);
-            document.getElementById("Awp").innerHTML = Awp;
+            document.getElementById("Awp").innerHTML = Awp.toFixed(3);
             document.getElementById("Cb").innerHTML = Cb.toFixed(3);
             document.getElementById("LCBx").innerHTML = LCBx.toFixed(3);
             document.getElementById("LCFx").innerHTML = LCFx.toFixed(3);
@@ -351,6 +387,13 @@ 
    Comparation
    
             document.getElementById("KG").innerHTML = KG.toFixed(3);
             document.getElementById("BM").innerHTML = BM.toFixed(3);
             document.getElementById("GM").innerHTML = GM.toFixed(3);
+            document.getElementById("trimd").innerHTML = trimd.toFixed(3);
+            document.getElementById("trima").innerHTML = trima.toFixed(3);
+            document.getElementById("heel").innerHTML = heel.toFixed(3);
+            document.getElementById("maxmom").innerHTML = maxmom.toFixed(3);
+            document.getElementById("shear").innerHTML = shear.toFixed(3);
+            document.getElementById("lightweight").innerHTML = lightweight.toFixed(3);
+            document.getElementById("deadweight").innerHTML = deadweight.toFixed(3);
             
             //calculated values
             let KBc = a.KB.toFixed(3);
@@ -359,13 +402,15 @@ 
    Comparation
    
             let draftc = ship.designState.calculationParameters.Draft_design;
             
             //baseObjects not sure why is not working
-            let lightweight = ship.structure.hull.attributes.lightweight + 15000;
-            let deadweight = 14000 + 56000;
+            let lightweightc = ship.structure.hull.attributes.lightweight + 15000;
+            let deadweightc = 14000 + 56000;
+            
             
             let lwlc = ship.structure.hull.calculateAttributesAtDraft(draftc).LWL;
             let bwlc = ship.structure.hull.calculateAttributesAtDraft(draftc).BWL;
             let Awpc = ship.structure.hull.calculateAttributesAtDraft(draftc).Awp;
             let cbc = ship.structure.hull.calculateAttributesAtDraft(draftc).Cb;
+            let dispc = lwlc * bwlc * cbc * draftc;
             
             document.getElementById("KBc").innerHTML = KBc;
             document.getElementById("KGc").innerHTML = KGc;
@@ -375,8 +420,9 @@ 
    Comparation
    
             document.getElementById("bwlc").innerHTML = bwlc.toFixed(3);
             document.getElementById("Awpc").innerHTML = Awpc.toFixed(3);
             document.getElementById("cbc").innerHTML = cbc.toFixed(3);
-            document.getElementById("lightweight").innerHTML = lightweight;
-            document.getElementById("deadweight").innerHTML = deadweight;
+            document.getElementById("lightweightc").innerHTML = lightweightc.toFixed(3);
+            document.getElementById("deadweightc").innerHTML = deadweightc.toFixed(3);
+            document.getElementById("dispc").innerHTML = dispc.toFixed(3);
             
             //% diff
             let KBd = (1-KB/KBc)*100;

From 043900a4e8dfc85d3bc19798ed4f9d63b529e86f Mon Sep 17 00:00:00 2001
From: EliasHasle 
Date: Fri, 3 Nov 2017 22:07:22 +0100
Subject: [PATCH 14/20] Simplification and correction of hull calculations, and
 more

Now it almost passes the blockCase_compare_json.html test. Only the z
component of the center of submerged volume deviates signfificantly from
the hand-calculated value. It looks like more zero-offset areas are sent
to the patchColumnCalculation function than I would expect. (Potential
bug.)

I have (temporarily) introduced simpler volume and area calculations
during debugging, and will try to reintroduce the bilinear patches, as
they give higher accuracy (at least for the volume).
---
 build/Vessel.js                        | 260 ++++++++++++-------------
 examples/blockCase_compare_json.html   |  29 +--
 source/CoreClasses/Hull.js             | 234 ++++++++++------------
 source/CoreClasses/Structure.js        |   2 +-
 source/functions/areaCalculations.js   |  14 +-
 source/functions/vectorOperations.js   |  11 ++
 source/functions/volumeCalculations.js |   9 +-
 7 files changed, 271 insertions(+), 288 deletions(-)

diff --git a/build/Vessel.js b/build/Vessel.js
index e638748..1f928c5 100644
--- a/build/Vessel.js
+++ b/build/Vessel.js
@@ -1,4 +1,4 @@
-//Vessel.js library, built 2017-10-27 15:17:23.678874, Checksum: d15ff7f6a612a6565c1b1e1ff341031a
+//Vessel.js library, built 2017-11-03 21:58:34.624540, Checksum: 66754104a0f843a0cff421b5aa24f15e
 /*
 Import like this in HTML:
 
@@ -33,11 +33,16 @@ function vecNormSquared(v) {
 	return v.x**2+v.y**2+v.z**2;
 }
 
+/*Adds two or more vectors given as individual parameters,
+and returns a new vector that is the component-wise 
+sum of the input vectors.*/
 function addVec(u,v, ...rest) {
 	if (rest.length > 0) return sumVec([u,v]+rest);
 	return {x: u.x+v.x, y: u.y+v.y, z: u.z+v.z};
 }
 
+//Takes an array of vectors as input, and returns a new vector
+//that is the component-wise sum of the input vectors.
 function sumVec(vectors) {
 	let S = {x:0, y:0, z:0};
 	for (let i = 0; i < vectors.length; i++) {
@@ -49,6 +54,12 @@ function sumVec(vectors) {
 	return S;
 }
 
+//Takes two vector parameters u,v, and returns the vector u-v.
+function subVec(u,v) {
+	//return addVec(u, scaleVec(v, -1)); //equivalent
+	return {x: u.x-v.x, y: u.y-v.y, z: u.z-v.z};
+}
+
 function dotProduct(u,v) {
 	return u.x*v.x + u.y*v.y + u.z*v.z;
 }
@@ -227,6 +238,7 @@ function combineAreas(array) {
 		yc /= A;
 	} else {
 		console.warn("Zero area combination.");
+		console.trace();
 		xc /= L;
 		yc /= L;
 	}
@@ -323,15 +335,20 @@ function bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22, segs=10) {
 	return A;
 }
 
-//Calculates the (arithmetic) average of the area of the two possible triangulations of the quad element (using two triangles).
+/*Calculates the (arithmetic) average of the area of the two possible triangulations of the quad element (using two triangles).
+This requires the base of the quad to be convex. If the base is arrowhead shaped,
+The calculation will fail in undefined ways.
+*/
 function elementArea(v1,v2,v3,v4) {
-	let A = 0.5*(Math.abs(signedTriangleArea(v1,v2,v3)) + Math.abs(signedTriangleArea(v2,v3,v4)) + Math.abs(signedTriangleArea(v3,v4,v1)) + Math.abs(signedTriangleArea(v4,v1,v2)));
+	let A1 = Math.abs(signedTriangleArea(v1,v2,v3)) + Math.abs(signedTriangleArea(v3,v4,v1));
+	let A2 = Math.abs(signedTriangleArea(v2,v3,v4)) + Math.abs(signedTriangleArea(v4,v1,v2));
+	let A = 0.5*(A1+A2);
 	return A;
 }
 
 function signedTriangleArea(v1,v2,v3) {
-	let u = addVec(v2,scaleVec(v1,-1));
-	let v = addVec(v3,scaleVec(v1,-1));
+	let u = subVec(v2,v1);
+	let v = subVec(v3,v1);
 	let c = crossProduct(u,v);
 	let A = 0.5*vecNorm(c);
 	return A;
@@ -342,18 +359,19 @@ function signedTriangleArea(v1,v2,v3) {
 //even though all x and y values are positive.
 //I am currently doing some tests here of an (over-)simplified calculation
 //The results so far indicate that, for the prism hull, the results are almost identical, except that with the simple calculation the center of volume is almost right (but wrong enough to disqualify such a simple calculation). The bug causing very wrong (too big) submerged volume is likely to be found elsewhere.
+/*Note that the coordinate system used here has xy as a grid, with z as heights on the grid, but in the intended application, which is calculations on transverse hull offsets, z corresponds to the vessel y axis, and y corresponds to the vessel z axis. In the application, the conversion between coordinate systems must be taken care of appropriately.*/
 											  // xy
 function patchColumnCalculation(x1, x2, y1, y2, z11, z12, z21, z22) {
 	//DEBUG START:
 	//Simpler approximate calculation of volume:
-	let z = 0.25*(z11+z12+z21+z22);
-	let V = Math.abs((x2-x1)*(y2-y1)*z);
+	let z = 0.25*(z11+z12+z21+z22); //"average height"
+	let V = Math.abs((x2-x1)*(y2-y1)*z); //base area times "average height"
 
 	//Very approximate center of volume
 	//(does not account for different signs on z values,
 	//but that should be OK for hull offsets)
-	let xc = (x1*(z11+z12)+x2*(z21+z22))/(z11+z12+z21+z22);
-	let yc = (y1*(z11+z21)+y2*(z12+z22))/(z11+z12+z21+z22);
+	let xc = (x1*(z11+z12)+x2*(z21+z22))/((z11+z12+z21+z22) || 1);
+	let yc = (y1*(z11+z21)+y2*(z12+z22))/((z11+z12+z21+z22) || 1);
 	let zc = 0.5*z;
 	
 	//Simple triangle average approximation for area
@@ -895,107 +913,49 @@ Object.assign(Hull.prototype, {
 		return output;
 	},
 	/*
-	This method is a mess, and introducing three different NaN correction has made it even worse in some ways. None of the correction modes fit for the usual cases. Maybe we should have a default mode that extrapolates to top NaNs and sets all other NaNs to zero. UPDATE: Now this is implemented as mode 3. Maybe we can skip the other ones.
+	Testing new version without nanCorrectionMode parameter, that defaults to setting lower NaNs to 0 and extrapolating highest data entry for upper NaNs (if existant, else set to 0). Inner NaNs will also be set to zero.
+	
+	This does not exactly work perfectly yet. getwaterline(0,3) gives interpolated values, even though the keel level is defined.
 	
 	Input:
 	z: level from bottom of ship (absolute value in meters)
-	nanCorrectionMode: 0 to set all NaNs to zero, 1 to output NaNs, 2 to replace NaNs with interpolated or extrapolated values, 3 to extrapolate to top NaNs and set all other NaNs to zero.
+	
+	Output:
+	Array representing waterline offsets for a given height from the keel (typically a draft).
 	*/
-	getWaterline: function(z, nanCorrectionMode=1) {
+	getWaterline: function(z) {
 		let ha = this.attributes;
-		let zr = z/ha.Depth;
-		let wls = this.halfBreadths.waterlines;
+		let zr = z/ha.Depth; //using zr requires fewer operations and less memory than a scaled copy of wls.
+		let wls = this.halfBreadths.waterlines;//.map(wl=>wl*ha.Depth);
 		let sts = this.halfBreadths.stations;
 		let tab = this.halfBreadths.table;
 
-		let {index: a, mu: mu} = bisectionSearch(wls, zr);
-		let wl;
-		if (a<0) {
-			if (nanCorrectionMode===0 || nanCorrectionMode===3) {
-				console.warn("getWaterLine: z below lowest defined waterline. Defaulting to zeros.");
-				return new Array(sts.length).fill(0);
-			}
-			if (nanCorrectionMode===1) {
-				console.warn("getWaterLine: z below lowest defined waterline. Outputting NaNs.");
-				return new Array(sts.length).fill(null);
-			}
-			else /*nanCorrectionMode===2*/ {
-				console.warn("getWaterLine: z below lowest defined waterline. Extrapolating lowest data entry.");
-				a=0;
-				mu=0;
-				//wl = tab[a].slice();
-			}			
-		} else if (a >= wls.length-1) {
-			if (nanCorrectionMode===0) {
-				console.warn("getWaterLine: z above highest defined waterline. Defaulting to zeros.");
+		if (zrwls[wls.length-1]) {
+				console.warn("getWaterLine: z above highest defined waterline. Proceeding with highest data entries.");
 				a = wls.length-2; //if this level is defined...
 				mu=1;
 				//wl = tab[a].slice();
-			}
-		}
-
-		//Linear interpolation between data waterlines
-		wl = new Array(sts.length);
-		for (let j = 0; j < wl.length; j++) {
-			if (nanCorrectionMode === 0) {
-				if (a+1 > wls.length-1) {
-					wl[j] = lerp(tab[a][j], 0, 0.5); //Suspicious!
-				} else {
-					wl[j] = lerp(tab[a][j] || 0, tab[a+1][j] || 0, mu || 0.5);
-				}
-			} else if (nanCorrectionMode === 1) {
-				if (a+1 > wls.length-1) {
-					wl[j] = lerp(tab[a][j], null, mu);
-				} else {
-					wl[j] = lerp(tab[a][j], tab[a+1][j], mu);
-				}
-			} else if (nanCorrectionMode === 2) {
-				//If necessary, sample from below
-				let b = a;
-				while (b>0 && isNaN(tab[b][j])) {
-					b--;
-				}
-				let lower;
-				if (b===0 && isNaN(tab[b][j])) {
-					lower = 0;
-				} else {
-					lower = tab[b][j];
-				}
-				//If necesary, sample from above
-				let c = a+1;
-				let upper;
-				if (c>wls.length-1) {
-					c = b;
-					upper = lower;
-				} else {
-					while (cwls.length-1 before the loop.
-					if (c===wls.length-1 && isNaN(tab[c][j])) {
-						//Fall back all the way to b
-						c = b;
-						upper = lower;
-					} else {
-						upper = tab[c][j];
-					}
+			} else {
+				({index: a, mu: mu} = bisectionSearch(wls, zr));
+				if (a === wls.length-1) {
+					a = wls.length-2;
+					mu = 1;
 				}
-				mu = c===b ? 0 : (a+(mu||0.5)-b)/(c-b); //what is this? Some kind of self-explanatory code, perhaps?
-				wl[j] = lerp(lower, upper, mu);
-			} else /*nanCorrectionMode === 3*/ {
+			}
+			
+			//Try to do linear interpolation between closest data waterlines, but handle null values well:
+			let wl = new Array(sts.length);
+			for (let j = 0; j < wl.length; j++) {
 				let lower, upper;
-				let b;
-				if (!isNaN(tab[a][j])) {
-					lower = tab[a][j];					
+				let b = a;
+				//Find lower value for interpolation
+				if (!isNaN(tab[b][j])) {
+					lower = tab[b][j];					
 				} else {
 					b = a+1;
 					while(b < wls.length && isNaN(tab[b][j])) {
@@ -1020,9 +980,11 @@ Object.assign(Hull.prototype, {
 						}
 					}
 				}
-				if (upper !== undefined) {}
-				else if (!isNaN(tab[a+1][j])) {
-					upper = tab[a+1][j];
+				//Find upper value for interpolation
+				let c = a+1;
+				if (upper !== undefined) {/*upper found above*/}
+				else if (!isNaN(tab[c][j])) {
+					upper = tab[c][j];
 				} else {
 					//The cell value is NaN.
 					//Upper is not defined.
@@ -1030,22 +992,22 @@ Object.assign(Hull.prototype, {
 					//or tab[a][j] is an inner NaN and
 					//there exists at least one number above it.
 					//In both cases I have to check above a+1.
-					b = a+2;
-					while (b < wls.length && isNaN(tab[b][j])) {
-						b++;
+					c = a+2;
+					while (c < wls.length && isNaN(tab[c][j])) {
+						c++;
 					}
-					if (b === wls.length) upper = lower;
+					if (c === wls.length) upper = lower;
 					else {
-						upper = tab[b][j];
+						upper = tab[c][j];
 					}
 				}
-				wl[j] = lerp(lower, upper, mu || 0.5);
+				//Linear interpolation
+				wl[j] = lerp(lower, upper, mu);
+						//Scale numerical values
+				if (!isNaN(wl[j])) wl[j] *= 0.5*ha.BOA;
 			}
-			
-			//Scale numerical values
-			if (!isNaN(wl[j])) wl[j] *= 0.5*ha.BOA;
-		}
 		return wl;
+		}
 	},
 	getStation: function(x) {
 		let ha = this.attributes;
@@ -1140,7 +1102,7 @@ Object.assign(Hull.prototype, {
 			}
 		}
 
-		//This does not yet account for undefined minY, maxY. Or does it?
+		//This does not yet account properly for undefined minY, maxY.
 		let port = [], star = [];
 		for (let i=0; ithis.attributes.Depth*wl);
 		let port = this.getStation(x);
+		if (!isNaN(maxZ)) {
+			let {index, mu} = bisectionSearch(wls, maxZ);
+			if (index < wls.length-1) {
+				wls[index+1] = lerp(wls[index], wls[index+1], mu);
+				port[index+1] = lerp(port[index], port[index+1], mu);
+				wls = wls.slice(0,index+2);
+				port = port.slice(0,index+2);
+			}
+		}
 		let star = port.map(hb=>-hb);
+
 		let sc = sectionCalculation({xs: wls, ymins: star, ymaxs: port});
 		return {
 			x: x, //or xc? or cg.. Hm.
@@ -1197,14 +1170,14 @@ Object.assign(Hull.prototype, {
 			A: sc.A,
 			Iz: sc.Ix,
 			Iy: sc.Iy,
-			maxX: sc.maxX,
-			minX: sc.minX,
+			maxZ: sc.maxX,
+			minZ: sc.minX,
 			maxY: sc.maxY,
 			minY: sc.minY
 		};
 	},
-	//Unoptimized, some redundant repetitions of calculations.
-	//NOT DONE YET. Outputs lots of NaN values.
+
+	//NOT DONE YET. Calculates too big Ap, Vs, Cb, and too small As for the test case. For the test, the Ap is exactly BWL*1m larger than it should be (why?). Too high Cb is clearly caused by too big Vs, and big Ap and Vs may have a common cause. The bilinear volume and area calculations have been temporarily replaced with simpler calculations, but this does not seem to help. I expect to find the bug(s) elsewhere.
 	//Important: calculateAttributesAtDraft takes one mandatory parameter T. (The function defined here is immediately called during construction of the prototype, and returns the proper function.)
 	calculateAttributesAtDraft: function() {
 		function levelCalculation(hull,
@@ -1228,22 +1201,43 @@ Object.assign(Hull.prototype, {
 			let lev = {};
 			Object.assign(lev, wlc);
 			//Projected area calculation (approximate):
-			lev.prMinY = wlc.minY || 0;
-			lev.prMaxY = wlc.maxY || 0;
-			lev.Ap = prev.Ap
-				+ trapezoidCalculation(prev.prMinY, prev.prMaxY, lev.prMinY, lev.prMaxY, prev.z, lev.z)["A"];
+			lev.prMinY = wlc.minY;
+			lev.prMaxY = wlc.maxY;
+			//DEBUG:
+			//console.info("prev.Ap = ", prev.Ap);
+			//console.info("Parameters to trapezoidCalculation: (%.2f, %.2f, %.2f, %.2f, %.2f, %.2f)", prev.prMinY, prev.prMaxY, lev.prMinY, lev.prMaxY, prev.z, z);
+			let AT = trapezoidCalculation(prev.prMinY, prev.prMaxY, lev.prMinY, lev.prMaxY, prev.z, z)["A"];
+			//console.log("Calculated area of trapezoid: ", AT);
+			lev.Ap = prev.Ap + AT;
+			//lev.Ap = prev.Ap
+			//	+ trapezoidCalculation(prev.prMinY, prev.prMaxY, lev.prMinY, lev.prMaxY, prev.z, z)["A"];
+			//DEBUG END
+
 			
 			//level bounds are for the bounding box of the submerged part of the hull
-			if (!isNaN(prev.minX) && prev.minX<=wlc.minX) 
+			if (!isNaN(wlc.minX) && wlc.minX<=prev.minX) 
+				lev.minX = wlc.minX;
+			else
 				lev.minX = prev.minX;
-			if (!isNaN(prev.maxX) && prev.maxX>=wlc.maxX) 
+			if (!isNaN(wlc.maxX) && wlc.maxX>=prev.maxX) 
+				lev.maxX = wlc.maxX;
+			else
 				lev.maxX = prev.maxX;
-			if (!isNaN(prev.minY) && prev.minY<=wlc.minY) 
+			if (!isNaN(wlc.minY) && wlc.minY<=prev.minY) 
+				lev.minY = wlc.minY;
+			else
 				lev.minY = prev.minY;
-			if (!isNaN(prev.maxY) && prev.maxY>=wlc.maxY) 
+			if (!isNaN(wlc.maxY) && wlc.maxY>=prev.maxY) 
+				lev.maxY = wlc.maxY;
+			else
 				lev.maxY = prev.maxY;
+			
 			lev.Vbb = (lev.maxX-lev.minX)*(lev.maxY-lev.minY)*z;
 			
+			//Keep level maxX and minX for finding end cap areas:
+			lev.maxXwp = wlc.maxX;
+			lev.minXwp = wlc.minX;
+			
 			//Find bilinear patches in the slice, and combine them.
 			//Many possibilities for getting the coordinate systems wrong.
 			let calculations = [];
@@ -1258,6 +1252,7 @@ Object.assign(Hull.prototype, {
 					patchColumnCalculation(sts[j], sts[j+1], prev.z, z, prwl[j], wl[j], prwl[j+1], wl[j+1]);
 				calculations.push(star);
 			}
+			console.log(calculations); //DEBUG
 			let C = combineVolumes(calculations);
 			//Cv of slice. Note that switching of yz must
 			//be done before combining with previous level
@@ -1266,8 +1261,11 @@ Object.assign(Hull.prototype, {
 			lev.Vs = prev.Vs + C.V; //hull volume below z
 			lev.As = prev.As + C.As; //outside surface below z
 
-			//Simple approximation (wrong) for end caps:
-			lev.As += 2*lev.Ap;
+			//End caps:
+			if (lev.minXwp <= sts[0])
+				lev.As += hull.stationCalculation(lev.minXwp, z)["A"];
+			if (lev.maxXwp >= sts[sts.length-1])
+				lev.As += hull.stationCalculation(lev.maxXwp, z)["A"];
 			
 			//center of volume below z (some potential for accumulated rounding error when calculating an accumulated average like this):
 			lev.Cv = scaleVec(addVec(
@@ -1302,12 +1300,14 @@ Object.assign(Hull.prototype, {
 				this.levelsNeedUpdate = false;
 			}
 			
-			//Find highest data waterline below water:
-			let {index: previ} = bisectionSearch(wls, T);
-			
-			//console.info("Highest data waterline below water: " + previ);
+			//Find highest data waterline below or at water level:
+			let {index, mu} = bisectionSearch(wls, T);
 			
-			let lc = levelCalculation(this, T, this.levels[previ] || undefined);
+			console.info("Highest data waterline below or at water level: " + index);
+			console.log(this.levels);
+			let lc;
+			if (mu===0) lc = this.levels[index];
+			else lc = levelCalculation(this, T, this.levels[index]);
 			
 			//Filter and rename for output
 			return {
diff --git a/examples/blockCase_compare_json.html b/examples/blockCase_compare_json.html
index 96dd71a..8021632 100644
--- a/examples/blockCase_compare_json.html
+++ b/examples/blockCase_compare_json.html
@@ -64,26 +64,6 @@
 	//logToScreen("Design state: ", ship.designState);
 
 	renderjson.set_show_to_level("all");
-		
-	//Hand-calculated values from @MrEranwe
-	//Some of these values are not calculated by the library yet.
-	//And the library does not account for trim when calculating Cb etc.
-	/*let disp = 102.4;
-	let draftms = 1.989;
-	//let draftap = 1.98;
-	//let draftfp = 1.999;
-	let lwl = 20;
-	let bwl = 5.046;
-	let Awet = 137.989;
-	let Awp = 99.95;
-	let Cb = disp / (draftms * lwl * bwl * 1.025);
-	let LCBx = 10.065;
-	let LCFx = 10.032;
-	let KB = 1.333;
-	let KG = 3.513;
-	let BM = 2.082;
-	let GM = KB + BM - KG;
-	//let trim = -0.019;*/
 	
 	let w = ship.getWeight();
 	logToScreen("Assuming weight from library calculation: ", w);
@@ -98,6 +78,11 @@
 	}, w);*/
 
 	
+	//Some of the hand-calculated values from @MrEranwe 
+	//are not calculated by the library yet.
+	//And the library does not account for trim when calculating Cb etc.
+	//Therefore I have made new, simpler, calculations that map to the ones made by the library.
+	
 	//Simplified calculations for triangular prism case, without trim:
 	let ha = ship.structure.hull.attributes;
 	let LOA = ha.LOA;
@@ -130,8 +115,8 @@
 	let maxYs = 0.5*BWL;
 	let LBP = LWL-ha.APP;
 	
-	let Tc = ship.calculateDraft();
-	logToScreen("Calculated draft: Man: "+ T + ", Machine: "+ Tc);
+	/*let Tc = ship.calculateDraft();
+	logToScreen("Calculated draft: Man: "+ T + ", Machine: "+ Tc);*/
 	logToScreen("Assuming draft = " + T + ":");
 	let hca = v.structure.hull.calculateAttributesAtDraft(T);
 	compToScreen("Calculated hull values at calculated draft: ", {
diff --git a/source/CoreClasses/Hull.js b/source/CoreClasses/Hull.js
index a6ab1a5..7505223 100644
--- a/source/CoreClasses/Hull.js
+++ b/source/CoreClasses/Hull.js
@@ -49,107 +49,47 @@ Object.assign(Hull.prototype, {
 		return output;
 	},
 	/*
-	This method is a mess, and introducing three different NaN correction has made it even worse in some ways. None of the correction modes fit for the usual cases. Maybe we should have a default mode that extrapolates to top NaNs and sets all other NaNs to zero. UPDATE: Now this is implemented as mode 3. Maybe we can skip the other ones.
+	Testing new version without nanCorrectionMode parameter, that defaults to setting lower NaNs to 0 and extrapolating highest data entry for upper NaNs (if existant, else set to 0). Inner NaNs will also be set to zero.
 	
 	Input:
 	z: level from bottom of ship (absolute value in meters)
-	nanCorrectionMode: 0 to set all NaNs to zero, 1 to output NaNs, 2 to replace NaNs with interpolated or extrapolated values, 3 to extrapolate to top NaNs and set all other NaNs to zero.
+	
+	Output:
+	Array representing waterline offsets for a given height from the keel (typically a draft).
 	*/
-	getWaterline: function(z, nanCorrectionMode=1) {
+	getWaterline: function(z) {
 		let ha = this.attributes;
-		let zr = z/ha.Depth;
-		let wls = this.halfBreadths.waterlines;
+		let zr = z/ha.Depth; //using zr requires fewer operations and less memory than a scaled copy of wls.
+		let wls = this.halfBreadths.waterlines;//.map(wl=>wl*ha.Depth);
 		let sts = this.halfBreadths.stations;
 		let tab = this.halfBreadths.table;
 
-		let {index: a, mu: mu} = bisectionSearch(wls, zr);
-		let wl;
-		if (a<0) {
-			if (nanCorrectionMode===0 || nanCorrectionMode===3) {
-				console.warn("getWaterLine: z below lowest defined waterline. Defaulting to zeros.");
+		if (zr= wls.length-1) {
-			if (nanCorrectionMode===0) {
-				console.warn("getWaterLine: z above highest defined waterline. Defaulting to zeros.");
-				return new Array(sts.length).fill(0);
-			}
-			if (nanCorrectionMode===1) {
-				console.warn("getWaterLine: z above highest defined waterline. Outputting NaNs.");
-				return new Array(sts.length).fill(null);
-			}
-			else /*nanCorrectionMode===2 || nanCorrectionMode===3*/ {
-				console.warn("getWaterLine: z above highest defined waterline. Proceeding with highest data entry.");
+		} else {
+			let a, mu;
+			if (zr>wls[wls.length-1]) {
+				console.warn("getWaterLine: z above highest defined waterline. Proceeding with highest data entries.");
 				a = wls.length-2; //if this level is defined...
 				mu=1;
 				//wl = tab[a].slice();
-			}
-		}
-
-		//Linear interpolation between data waterlines
-		wl = new Array(sts.length);
-		for (let j = 0; j < wl.length; j++) {
-			if (nanCorrectionMode === 0) {
-				if (a+1 > wls.length-1) {
-					wl[j] = lerp(tab[a][j], 0, 0.5); //Suspicious!
-				} else {
-					wl[j] = lerp(tab[a][j] || 0, tab[a+1][j] || 0, mu || 0.5);
-				}
-			} else if (nanCorrectionMode === 1) {
-				if (a+1 > wls.length-1) {
-					wl[j] = lerp(tab[a][j], null, mu);
-				} else {
-					wl[j] = lerp(tab[a][j], tab[a+1][j], mu);
-				}
-			} else if (nanCorrectionMode === 2) {
-				//If necessary, sample from below
-				let b = a;
-				while (b>0 && isNaN(tab[b][j])) {
-					b--;
-				}
-				let lower;
-				if (b===0 && isNaN(tab[b][j])) {
-					lower = 0;
-				} else {
-					lower = tab[b][j];
-				}
-				//If necesary, sample from above
-				let c = a+1;
-				let upper;
-				if (c>wls.length-1) {
-					c = b;
-					upper = lower;
-				} else {
-					while (cwls.length-1 before the loop.
-					if (c===wls.length-1 && isNaN(tab[c][j])) {
-						//Fall back all the way to b
-						c = b;
-						upper = lower;
-					} else {
-						upper = tab[c][j];
-					}
+			} else {
+				({index: a, mu: mu} = bisectionSearch(wls, zr));
+				if (a === wls.length-1) {
+					a = wls.length-2;
+					mu = 1;
 				}
-				mu = c===b ? 0 : (a+(mu||0.5)-b)/(c-b); //what is this? Some kind of self-explanatory code, perhaps?
-				wl[j] = lerp(lower, upper, mu);
-			} else /*nanCorrectionMode === 3*/ {
+			}
+			
+			//Try to do linear interpolation between closest data waterlines, but handle null values well:
+			let wl = new Array(sts.length);
+			for (let j = 0; j < wl.length; j++) {
 				let lower, upper;
-				let b;
-				if (!isNaN(tab[a][j])) {
-					lower = tab[a][j];					
+				let b = a;
+				//Find lower value for interpolation
+				if (!isNaN(tab[b][j])) {
+					lower = tab[b][j];					
 				} else {
 					b = a+1;
 					while(b < wls.length && isNaN(tab[b][j])) {
@@ -174,9 +114,11 @@ Object.assign(Hull.prototype, {
 						}
 					}
 				}
-				if (upper !== undefined) {}
-				else if (!isNaN(tab[a+1][j])) {
-					upper = tab[a+1][j];
+				//Find upper value for interpolation
+				let c = a+1;
+				if (upper !== undefined) {/*upper found above*/}
+				else if (!isNaN(tab[c][j])) {
+					upper = tab[c][j];
 				} else {
 					//The cell value is NaN.
 					//Upper is not defined.
@@ -184,22 +126,22 @@ Object.assign(Hull.prototype, {
 					//or tab[a][j] is an inner NaN and
 					//there exists at least one number above it.
 					//In both cases I have to check above a+1.
-					b = a+2;
-					while (b < wls.length && isNaN(tab[b][j])) {
-						b++;
+					c = a+2;
+					while (c < wls.length && isNaN(tab[c][j])) {
+						c++;
 					}
-					if (b === wls.length) upper = lower;
+					if (c === wls.length) upper = lower;
 					else {
-						upper = tab[b][j];
+						upper = tab[c][j];
 					}
 				}
-				wl[j] = lerp(lower, upper, mu || 0.5);
+				//Linear interpolation
+				wl[j] = lerp(lower, upper, mu);
+						//Scale numerical values
+				if (!isNaN(wl[j])) wl[j] *= 0.5*ha.BOA;
 			}
-			
-			//Scale numerical values
-			if (!isNaN(wl[j])) wl[j] *= 0.5*ha.BOA;
-		}
 		return wl;
+		}
 	},
 	getStation: function(x) {
 		let ha = this.attributes;
@@ -235,13 +177,13 @@ Object.assign(Hull.prototype, {
 	//THIS is a candidate for causing wrong Ix, Iy values.
 	//Much logic that can go wrong.
 									//typically deck bounds
-	waterlineCalculation: function(z, bounds, nanCorrectionMode=3) {
+	waterlineCalculation: function(z, bounds) {
 		let {minX, maxX, minY, maxY} = bounds || {};
 
 		console.group/*Collapsed*/("waterlineCalculation.");
-		console.info("Arguments: z=", z, " Boundaries: ", arguments[1], " NaN correction mode: ", nanCorrectionMode);
+		console.info("Arguments: z=", z, " Boundaries: ", arguments[1]);
 		
-		let wl = this.getWaterline(z, nanCorrectionMode);
+		let wl = this.getWaterline(z);
 		console.info("wl: ", wl); //DEBUG
 
 		let LOA = this.attributes.LOA;
@@ -294,7 +236,7 @@ Object.assign(Hull.prototype, {
 			}
 		}
 
-		//This does not yet account for undefined minY, maxY. Or does it?
+		//This does not yet account properly for undefined minY, maxY.
 		let port = [], star = [];
 		for (let i=0; ithis.attributes.Depth*wl);
 		let port = this.getStation(x);
+		if (!isNaN(maxZ)) {
+			let {index, mu} = bisectionSearch(wls, maxZ);
+			if (index < wls.length-1) {
+				wls[index+1] = lerp(wls[index], wls[index+1], mu);
+				port[index+1] = lerp(port[index], port[index+1], mu);
+				wls = wls.slice(0,index+2);
+				port = port.slice(0,index+2);
+			}
+		}
 		let star = port.map(hb=>-hb);
+
 		let sc = sectionCalculation({xs: wls, ymins: star, ymaxs: port});
 		return {
 			x: x, //or xc? or cg.. Hm.
@@ -351,14 +304,14 @@ Object.assign(Hull.prototype, {
 			A: sc.A,
 			Iz: sc.Ix,
 			Iy: sc.Iy,
-			maxX: sc.maxX,
-			minX: sc.minX,
+			maxZ: sc.maxX,
+			minZ: sc.minX,
 			maxY: sc.maxY,
 			minY: sc.minY
 		};
 	},
-	//Unoptimized, some redundant repetitions of calculations.
-	//NOT DONE YET. Outputs lots of NaN values.
+
+	//NOT DONE YET. Calculates too big Ap, Vs, Cb, and too small As for the test case. For the test, the Ap is exactly BWL*1m larger than it should be (why?). Too high Cb is clearly caused by too big Vs, and big Ap and Vs may have a common cause. The bilinear volume and area calculations have been temporarily replaced with simpler calculations, but this does not seem to help. I expect to find the bug(s) elsewhere.
 	//Important: calculateAttributesAtDraft takes one mandatory parameter T. (The function defined here is immediately called during construction of the prototype, and returns the proper function.)
 	calculateAttributesAtDraft: function() {
 		function levelCalculation(hull,
@@ -378,32 +331,53 @@ Object.assign(Hull.prototype, {
 				Cv: {x:0, y:0, z:0}
 			}) {
 			
-			let wlc = hull.waterlineCalculation(z,{},3);
+			let wlc = hull.waterlineCalculation(z,{});
 			let lev = {};
 			Object.assign(lev, wlc);
 			//Projected area calculation (approximate):
-			lev.prMinY = wlc.minY || 0;
-			lev.prMaxY = wlc.maxY || 0;
-			lev.Ap = prev.Ap
-				+ trapezoidCalculation(prev.prMinY, prev.prMaxY, lev.prMinY, lev.prMaxY, prev.z, lev.z)["A"];
+			lev.prMinY = wlc.minY;
+			lev.prMaxY = wlc.maxY;
+			//DEBUG:
+			//console.info("prev.Ap = ", prev.Ap);
+			//console.info("Parameters to trapezoidCalculation: (%.2f, %.2f, %.2f, %.2f, %.2f, %.2f)", prev.prMinY, prev.prMaxY, lev.prMinY, lev.prMaxY, prev.z, z);
+			let AT = trapezoidCalculation(prev.prMinY, prev.prMaxY, lev.prMinY, lev.prMaxY, prev.z, z)["A"];
+			//console.log("Calculated area of trapezoid: ", AT);
+			lev.Ap = prev.Ap + AT;
+			//lev.Ap = prev.Ap
+			//	+ trapezoidCalculation(prev.prMinY, prev.prMaxY, lev.prMinY, lev.prMaxY, prev.z, z)["A"];
+			//DEBUG END
+
 			
 			//level bounds are for the bounding box of the submerged part of the hull
-			if (!isNaN(prev.minX) && prev.minX<=wlc.minX) 
+			if (!isNaN(wlc.minX) && wlc.minX<=prev.minX) 
+				lev.minX = wlc.minX;
+			else
 				lev.minX = prev.minX;
-			if (!isNaN(prev.maxX) && prev.maxX>=wlc.maxX) 
+			if (!isNaN(wlc.maxX) && wlc.maxX>=prev.maxX) 
+				lev.maxX = wlc.maxX;
+			else
 				lev.maxX = prev.maxX;
-			if (!isNaN(prev.minY) && prev.minY<=wlc.minY) 
+			if (!isNaN(wlc.minY) && wlc.minY<=prev.minY) 
+				lev.minY = wlc.minY;
+			else
 				lev.minY = prev.minY;
-			if (!isNaN(prev.maxY) && prev.maxY>=wlc.maxY) 
+			if (!isNaN(wlc.maxY) && wlc.maxY>=prev.maxY) 
+				lev.maxY = wlc.maxY;
+			else
 				lev.maxY = prev.maxY;
+			
 			lev.Vbb = (lev.maxX-lev.minX)*(lev.maxY-lev.minY)*z;
 			
+			//Keep level maxX and minX for finding end cap areas:
+			lev.maxXwp = wlc.maxX;
+			lev.minXwp = wlc.minX;
+			
 			//Find bilinear patches in the slice, and combine them.
 			//Many possibilities for getting the coordinate systems wrong.
 			let calculations = [];
 			let sts = hull.halfBreadths.stations.map(st=>st*hull.attributes.LOA);
-			let wl = hull.getWaterline(z,3);
-			let prwl = hull.getWaterline(prev.z,3);
+			let wl = hull.getWaterline(z);
+			let prwl = hull.getWaterline(prev.z);
 			for (let j = 0; j < sts.length-1; j++) {
 				let port = 
 					patchColumnCalculation(sts[j], sts[j+1], prev.z, z, -prwl[j], -wl[j], -prwl[j+1], -wl[j+1]);
@@ -412,6 +386,7 @@ Object.assign(Hull.prototype, {
 					patchColumnCalculation(sts[j], sts[j+1], prev.z, z, prwl[j], wl[j], prwl[j+1], wl[j+1]);
 				calculations.push(star);
 			}
+			console.log(calculations); //DEBUG
 			let C = combineVolumes(calculations);
 			//Cv of slice. Note that switching of yz must
 			//be done before combining with previous level
@@ -420,8 +395,11 @@ Object.assign(Hull.prototype, {
 			lev.Vs = prev.Vs + C.V; //hull volume below z
 			lev.As = prev.As + C.As; //outside surface below z
 
-			//Simple (and wrong) approximation for end caps:
-			lev.As += 2*lev.Ap; //(will work for prism case)
+			//End caps:
+			if (lev.minXwp <= sts[0])
+				lev.As += hull.stationCalculation(lev.minXwp, z)["A"];
+			if (lev.maxXwp >= sts[sts.length-1])
+				lev.As += hull.stationCalculation(lev.maxXwp, z)["A"];
 			
 			//center of volume below z (some potential for accumulated rounding error when calculating an accumulated average like this):
 			lev.Cv = scaleVec(addVec(
@@ -456,12 +434,14 @@ Object.assign(Hull.prototype, {
 				this.levelsNeedUpdate = false;
 			}
 			
-			//Find highest data waterline below water:
-			let {index: previ} = bisectionSearch(wls, T);
-			
-			//console.info("Highest data waterline below water: " + previ);
+			//Find highest data waterline below or at water level:
+			let {index, mu} = bisectionSearch(wls, T);
 			
-			let lc = levelCalculation(this, T, this.levels[previ] || undefined);
+			console.info("Highest data waterline below or at water level: " + index);
+			console.log(this.levels);
+			let lc;
+			if (mu===0) lc = this.levels[index];
+			else lc = levelCalculation(this, T, this.levels[index]);
 			
 			//Filter and rename for output
 			return {
diff --git a/source/CoreClasses/Structure.js b/source/CoreClasses/Structure.js
index 26b677b..083c493 100644
--- a/source/CoreClasses/Structure.js
+++ b/source/CoreClasses/Structure.js
@@ -63,7 +63,7 @@ Object.assign(Structure.prototype, {
 			let zc = d.zFloor+0.5*d.thickness;
 			let yc = d.yCentre;
 			let b = d.breadth;
-			let wlc = this.hull.waterlineCalculation(zc, {minX: d.xAft, maxX: d.xFwd, minY: yc-0.5*b, maxY: yc+0.5*b}, 3);
+			let wlc = this.hull.waterlineCalculation(zc, {minX: d.xAft, maxX: d.xFwd, minY: yc-0.5*b, maxY: yc+0.5*b});
 			components.push({
 				//approximation
 				mass: wlc.Awp*d.thickness*d.density,
diff --git a/source/functions/areaCalculations.js b/source/functions/areaCalculations.js
index cda838e..0eb500a 100644
--- a/source/functions/areaCalculations.js
+++ b/source/functions/areaCalculations.js
@@ -70,6 +70,7 @@ function combineAreas(array) {
 		yc /= A;
 	} else {
 		console.warn("Zero area combination.");
+		console.trace();
 		xc /= L;
 		yc /= L;
 	}
@@ -166,15 +167,20 @@ function bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22, segs=10) {
 	return A;
 }
 
-//Calculates the (arithmetic) average of the area of the two possible triangulations of the quad element (using two triangles).
+/*Calculates the (arithmetic) average of the area of the two possible triangulations of the quad element (using two triangles).
+This requires the base of the quad to be convex. If the base is arrowhead shaped,
+The calculation will fail in undefined ways.
+*/
 function elementArea(v1,v2,v3,v4) {
-	let A = 0.5*(Math.abs(signedTriangleArea(v1,v2,v3)) + Math.abs(signedTriangleArea(v2,v3,v4)) + Math.abs(signedTriangleArea(v3,v4,v1)) + Math.abs(signedTriangleArea(v4,v1,v2)));
+	let A1 = Math.abs(signedTriangleArea(v1,v2,v3)) + Math.abs(signedTriangleArea(v3,v4,v1));
+	let A2 = Math.abs(signedTriangleArea(v2,v3,v4)) + Math.abs(signedTriangleArea(v4,v1,v2));
+	let A = 0.5*(A1+A2);
 	return A;
 }
 
 function signedTriangleArea(v1,v2,v3) {
-	let u = addVec(v2,scaleVec(v1,-1));
-	let v = addVec(v3,scaleVec(v1,-1));
+	let u = subVec(v2,v1);
+	let v = subVec(v3,v1);
 	let c = crossProduct(u,v);
 	let A = 0.5*vecNorm(c);
 	return A;
diff --git a/source/functions/vectorOperations.js b/source/functions/vectorOperations.js
index 27227f7..c6c32dc 100644
--- a/source/functions/vectorOperations.js
+++ b/source/functions/vectorOperations.js
@@ -21,11 +21,16 @@ function vecNormSquared(v) {
 	return v.x**2+v.y**2+v.z**2;
 }
 
+/*Adds two or more vectors given as individual parameters,
+and returns a new vector that is the component-wise 
+sum of the input vectors.*/
 function addVec(u,v, ...rest) {
 	if (rest.length > 0) return sumVec([u,v]+rest);
 	return {x: u.x+v.x, y: u.y+v.y, z: u.z+v.z};
 }
 
+//Takes an array of vectors as input, and returns a new vector
+//that is the component-wise sum of the input vectors.
 function sumVec(vectors) {
 	let S = {x:0, y:0, z:0};
 	for (let i = 0; i < vectors.length; i++) {
@@ -37,6 +42,12 @@ function sumVec(vectors) {
 	return S;
 }
 
+//Takes two vector parameters u,v, and returns the vector u-v.
+function subVec(u,v) {
+	//return addVec(u, scaleVec(v, -1)); //equivalent
+	return {x: u.x-v.x, y: u.y-v.y, z: u.z-v.z};
+}
+
 function dotProduct(u,v) {
 	return u.x*v.x + u.y*v.y + u.z*v.z;
 }
diff --git a/source/functions/volumeCalculations.js b/source/functions/volumeCalculations.js
index e036d33..1f8b321 100644
--- a/source/functions/volumeCalculations.js
+++ b/source/functions/volumeCalculations.js
@@ -5,18 +5,19 @@
 //even though all x and y values are positive.
 //I am currently doing some tests here of an (over-)simplified calculation
 //The results so far indicate that, for the prism hull, the results are almost identical, except that with the simple calculation the center of volume is almost right (but wrong enough to disqualify such a simple calculation). The bug causing very wrong (too big) submerged volume is likely to be found elsewhere.
+/*Note that the coordinate system used here has xy as a grid, with z as heights on the grid, but in the intended application, which is calculations on transverse hull offsets, z corresponds to the vessel y axis, and y corresponds to the vessel z axis. In the application, the conversion between coordinate systems must be taken care of appropriately.*/
 											  // xy
 function patchColumnCalculation(x1, x2, y1, y2, z11, z12, z21, z22) {
 	//DEBUG START:
 	//Simpler approximate calculation of volume:
-	let z = 0.25*(z11+z12+z21+z22);
-	let V = Math.abs((x2-x1)*(y2-y1)*z);
+	let z = 0.25*(z11+z12+z21+z22); //"average height"
+	let V = Math.abs((x2-x1)*(y2-y1)*z); //base area times "average height"
 
 	//Very approximate center of volume
 	//(does not account for different signs on z values,
 	//but that should be OK for hull offsets)
-	let xc = (x1*(z11+z12)+x2*(z21+z22))/(z11+z12+z21+z22);
-	let yc = (y1*(z11+z21)+y2*(z12+z22))/(z11+z12+z21+z22);
+	let xc = (x1*(z11+z12)+x2*(z21+z22))/((z11+z12+z21+z22) || 1);
+	let yc = (y1*(z11+z21)+y2*(z12+z22))/((z11+z12+z21+z22) || 1);
 	let zc = 0.5*z;
 	
 	//Simple triangle average approximation for area

From 63b0a7c51c32b2192d217fcb205903c25a57b58d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=8Dcaro=20Fonseca?= 
Date: Wed, 8 Nov 2017 16:42:35 +0100
Subject: [PATCH 15/20] Hull Hydrostatics example

Work in progress.
---
 examples/Hull_hydrostatics.html               | 288 ++++++++++++++++++
 .../ship_specifications/prysmaticHull.json    |  68 +++++
 2 files changed, 356 insertions(+)
 create mode 100644 examples/Hull_hydrostatics.html
 create mode 100644 examples/data/ship_specifications/prysmaticHull.json

diff --git a/examples/Hull_hydrostatics.html b/examples/Hull_hydrostatics.html
new file mode 100644
index 0000000..e9e7729
--- /dev/null
+++ b/examples/Hull_hydrostatics.html
@@ -0,0 +1,288 @@
+
+
+
+
+
+
+    
+    
+    
+    
+    
+
+    Hull Hydrostatics
+
+    
+    
+    
+    
+    
+
+    
+     
+	
+	
+	
+	
+	
+    
+	
+    
+
+	
+	
+
+	
+	
+	
+	
+	
+
+    
+    
+    
+
+
+
+
+
+    
+    
+            
+    
+    
    
+        
+        
+        
    
+            
+            
+            
+            
+             
    
+                
    Hull Hydrostatics
    
+            
      
+            
+			
    
+                
+                
+                
+                
    Interactive calculation of hull hydrostatic properties.
    
+				                                    
+				
    Developed by: Elias Hasle, Vicente Alejandro Iváñez Encinas and Henrique M. Gaspar. Visualization made using three.js.
    
+
+            
    
+			
+            
    
+                
    Input
    
+                
+                
    Contents of the file:
    
+                
    
+                
    
+            
    
+            
+            
    
+                
    3D orbit view of hull
    
+			
    
+			
    
+            
    
+                
+			
    
+            
    
+                
    Hydrostatic Properties
    
+                
+                    
+                        
+                            
+                            
+                        
+                    
+                    
+                    
+                        
+                            
+                            
+                        
+                        
+                            
+                            
+                        
+                        
+                           
+						    
+                        
+                        
+                           
+						    
+                        
+                        
+                           
+						   
+                        
+                        
+                           
+						    
+                        
+                        
+                           
+						    
+                        
+                        
+                           
+						    
+                        
+                        
+                        
+                    
+                
+                
    VariableResults
    Draft Amidships (m)
    Displacement (t)
    Waterline Length (m)
    Maximum Waterline Breadth (m)
    Waterplane Area (m2)
    Cb
    KB (m)
    BM (m)
    
+			
    
+        
    
+    
+        
+    
+    
+    
+    
+    
+    
+      
+
+
+
+
\ No newline at end of file
diff --git a/examples/data/ship_specifications/prysmaticHull.json b/examples/data/ship_specifications/prysmaticHull.json
new file mode 100644
index 0000000..8968b94
--- /dev/null
+++ b/examples/data/ship_specifications/prysmaticHull.json
@@ -0,0 +1,68 @@
+{
+  "attributes": {},
+  "designState": {
+    "calculationParameters": {
+      "LWL_design": 20,
+      "Draft_design": 2,
+      "Cb_design": 0.5,
+      "speed": 12,
+      "crew": 20,
+      "K": 0.032,
+      "MCR": 10000,
+      "SFR": 0.000215,
+      "Co": 0.3,
+      "tripDuration": 40
+    },
+    "objectOverrides": {
+      "common": {
+        "fullness": 0.7
+      }
+    }
+  },
+  "structure": {
+    "hull": {
+      "attributes": {
+        "LOA": 20,
+        "BOA": 10,
+        "Depth": 4,
+        "APP": 0,
+        "lightweight": 15000
+      },
+      "halfBreadths": {
+        "waterlines": [
+          0,
+          1
+        ],
+        "stations": [
+          0,
+          1
+        ],
+        "table": [
+          [
+            0,
+            0
+          ],
+          [
+            1,
+            1
+          ]
+        ]
+      },
+      "buttockHeights": {}
+    },
+    "decks": {
+      "WheelHouseTop": {
+        "zFloor": 4,
+        "thickness": 0.3,
+        "xAft": 0,
+        "xFwd": 20,
+        "yCentre": 0,
+        "breadth": 10,
+        "density": 1500
+      }
+    },
+    "bulkheads": {}
+  },
+  "baseObjects": [],
+  "derivedObjects": []
+}

From 53b4fe4fa035f322aaba299934e1bf4b617ed8c9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=8Dcaro=20Fonseca?= 
Date: Wed, 8 Nov 2017 18:19:27 +0100
Subject: [PATCH 16/20] Cotinuing example

---
 examples/Hull_hydrostatics.html | 51 +++++++++++++++++++++++----------
 1 file changed, 36 insertions(+), 15 deletions(-)

diff --git a/examples/Hull_hydrostatics.html b/examples/Hull_hydrostatics.html
index e9e7729..a0dd057 100644
--- a/examples/Hull_hydrostatics.html
+++ b/examples/Hull_hydrostatics.html
@@ -1,4 +1,5 @@
 
+
 
 
 
@@ -23,6 +24,9 @@
 	
 	
 	
+
+	
+	
 	
 	
     
@@ -161,6 +165,7 @@ 
    Hydrostatic Properties
    
 		
 	    // Create scene+
 		var renderer, scene, camera, controls, ship3D;
+		var gui = new dat.GUI();
 		
 		//Ready renderer and scene
         (function (){
@@ -234,23 +239,39 @@ 
    Hydrostatic Properties
    
 			camera.position.set(0.7*LOA, 0.7*LOA, 0.7*LOA);
 			controls.target = new THREE.Vector3(LOA/2,0,0);
 			controls.update();
-			animate();
-            
+
             //calculate hydrostatics
 			let {KB: KBc, BM: BMc} = ship.calculateStability();
-            let draftc = ship.designState.calculationParameters.Draft_design;
-            draftc = 1;
-			let {LWL: lwlc, Cb: cbc, BWL: bwlc, Awp: Awpc} = ship.structure.hull.calculateAttributesAtDraft(draftc);
-            let dispc = lwlc * bwlc * cbc * draftc;
-			
-            document.getElementById("KBc").innerHTML = KBc.toFixed(3);
-            document.getElementById("BMc").innerHTML = BMc.toFixed(3);
-            document.getElementById("draftc").innerHTML = draftc;
-            document.getElementById("lwlc").innerHTML = lwlc.toFixed(3);
-            document.getElementById("bwlc").innerHTML = bwlc.toFixed(3);
-            document.getElementById("Awpc").innerHTML = Awpc.toFixed(3);
-            document.getElementById("cbc").innerHTML = cbc.toFixed(3);
-            document.getElementById("dispc").innerHTML = dispc.toFixed(3);
+			let {LWL: lwlc, Cb: cbc, BWL: bwlc, Awp: Awpc} = ship.structure.hull.calculateAttributesAtDraft(ship.designState.calculationParameters.Draft_design);
+			let dispc = lwlc * bwlc * cbc * ship.designState.calculationParameters.Draft_design;
+
+			document.getElementById("KBc").innerHTML = KBc.toFixed(3);
+			document.getElementById("BMc").innerHTML = BMc.toFixed(3);
+			document.getElementById("draftc").innerHTML = ship.designState.calculationParameters.Draft_design.toFixed(3);
+			document.getElementById("lwlc").innerHTML = lwlc.toFixed(3);
+			document.getElementById("bwlc").innerHTML = bwlc.toFixed(3);
+			document.getElementById("Awpc").innerHTML = Awpc.toFixed(3);
+			document.getElementById("cbc").innerHTML = cbc.toFixed(3);
+			document.getElementById("dispc").innerHTML = dispc.toFixed(3);
+
+			gui.add(ship.designState.calculationParameters, "Draft_design", 0, 4, 0.01).onLoad(function (value) {
+				let {KB: KBc, BM: BMc} = ship.calculateStability();
+				let {LWL: lwlc, Cb: cbc, BWL: bwlc, Awp: Awpc} = ship.structure.hull.calculateAttributesAtDraft(ship.designState.calculationParameters.Draft_design);
+				let dispc = lwlc * bwlc * cbc * ship.designState.calculationParameters.Draft_design;
+
+				document.getElementById("KBc").innerHTML = KBc.toFixed(3);
+				document.getElementById("BMc").innerHTML = BMc.toFixed(3);
+				document.getElementById("draftc").innerHTML = ship.designState.calculationParameters.Draft_design.toFixed(3);
+				document.getElementById("lwlc").innerHTML = lwlc.toFixed(3);
+				document.getElementById("bwlc").innerHTML = bwlc.toFixed(3);
+				document.getElementById("Awpc").innerHTML = Awpc.toFixed(3);
+				document.getElementById("cbc").innerHTML = cbc.toFixed(3);
+				document.getElementById("dispc").innerHTML = dispc.toFixed(3);
+				ship3D.position.z = -ship.designState.calculationParameters.Draft_design;
+
+			});
+
+			animate();
 			
 		}
 		

From 949f958cc5125107706c1af9529e036902ee631f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=8Dcaro=20Fonseca?= 
Date: Wed, 8 Nov 2017 18:23:55 +0100
Subject: [PATCH 17/20] Minor correction

---
 examples/Hull_hydrostatics.html | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/Hull_hydrostatics.html b/examples/Hull_hydrostatics.html
index a0dd057..caebeeb 100644
--- a/examples/Hull_hydrostatics.html
+++ b/examples/Hull_hydrostatics.html
@@ -254,7 +254,7 @@ 
    Hydrostatic Properties
    
 			document.getElementById("cbc").innerHTML = cbc.toFixed(3);
 			document.getElementById("dispc").innerHTML = dispc.toFixed(3);
 
-			gui.add(ship.designState.calculationParameters, "Draft_design", 0, 4, 0.01).onLoad(function (value) {
+			gui.add(ship.designState.calculationParameters, "Draft_design", 0, 4, 0.01).onChange(function (value) {
 				let {KB: KBc, BM: BMc} = ship.calculateStability();
 				let {LWL: lwlc, Cb: cbc, BWL: bwlc, Awp: Awpc} = ship.structure.hull.calculateAttributesAtDraft(ship.designState.calculationParameters.Draft_design);
 				let dispc = lwlc * bwlc * cbc * ship.designState.calculationParameters.Draft_design;

From c0535c6f6eaeb3b282e7498933fcd1157689da28 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=8Dcaro=20Fonseca?= 
Date: Wed, 8 Nov 2017 18:53:40 +0100
Subject: [PATCH 18/20] Final touches to hydrostatics example

---
 examples/Hull_hydrostatics.html | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/examples/Hull_hydrostatics.html b/examples/Hull_hydrostatics.html
index caebeeb..e107c89 100644
--- a/examples/Hull_hydrostatics.html
+++ b/examples/Hull_hydrostatics.html
@@ -144,11 +144,11 @@ 
    Hydrostatic Properties
    
 						    
                         
                         
-                           KB (m)
+                           KB (m) WRONG VALUE
 						    
                         
                         
-                           BM (m)
+                           BM (m) WRONG VALUE
 						    
                         
                         
@@ -164,7 +164,7 @@ 
    Hydrostatic Properties
    
 		"use strict";
 		
 	    // Create scene+
-		var renderer, scene, camera, controls, ship3D;
+		var renderer, scene, camera, controls, ship3D, ocean;
 		var gui = new dat.GUI();
 		
 		//Ready renderer and scene
@@ -203,9 +203,17 @@ 
    Hydrostatic Properties
    
 				sun.position.set(1,1,1);
 				return sun;
 			}());
+
+			//Add rudimentary ocean, to see effect of draft changes:
+			ocean = new THREE.Mesh(new THREE.PlaneBufferGeometry(1000,1000), new THREE.MeshPhongMaterial({color: 0x2222bb}));
+			scene.add(ocean);
+			var grid = new THREE.GridHelper(1000, 100);
+			grid.rotation.x = 0.5*Math.PI;
+			grid.position.z = 0.01;
+			scene.add(grid);
 			
         })();
-				
+		
         // read file from user
         function readSingleFile(e){
             var file = e.target.files[0];

From e23c73a5f0c822cd73fce24cc753a588e0f848d7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Vicente=20Alejandro=20Iv=C3=A1=C3=B1ez=20Encinas?=
 
Date: Thu, 9 Nov 2017 19:19:13 +0100
Subject: [PATCH 19/20] Add different tables & "run around" for baseObject
 problem

---
 examples/blockCase_compare.html | 218 +++++++++++++++++++++++++++-----
 1 file changed, 184 insertions(+), 34 deletions(-)

diff --git a/examples/blockCase_compare.html b/examples/blockCase_compare.html
index a59685d..3549f23 100644
--- a/examples/blockCase_compare.html
+++ b/examples/blockCase_compare.html
@@ -105,12 +105,46 @@ 
    3D orbit view of hull and parts
    
             
             
             
    
-                
    Calculations
    
-                
    GM : m.
    
+                
    Whole ship
    
+                
    Ship specifications
    
+                
+                    
+                        
+                            
+                            
+                        
+                    
+                    
+                        
+                            
+                            
+                        
+                        
+                            
+                            
+                        
+                        
+                            
+                            
+                        
+                        
+                            
+                            
+                        
+                        
+                            
+                            
+                        
+                        
+                            
+                            
+                        
+                        
+                    
                 
-			
-            
    
-                
    Comparation
    
+
+                
    VariableValue
    Lenght Overall m
    Breadth Overall m
    Depth m
    Design Draft m
    Design speed knots
    Lightweight kg
    
+                
    Calculations
    
                 
@@ -144,7 +178,7 @@ 
    Comparation
     
-                            
+                            
@@ -174,7 +208,7 @@ 
    Comparation
     
-                            
+                            
@@ -211,18 +245,6 @@ 
    Comparation
     
-                        
-                         
-                            
-                            
-                            
-                        
-                        
-                            
-                            
-                            
@@ -258,6 +280,74 @@ 
    Comparation
    
                     
                         
    
                            Displacement (t) 
                             
                              
                         
    
                         
    
                            Waterline Lenght (m)
                            Cb 
                             
                             
                         
    
                         
    
                            LCBx (m)
                             
                            
                            
    Lightweigth (kg)
-                           
    Deadweight (kg)
-                            
                         
    
                         
    
                            Trim (m)
                     
                 
                 
    
+                
    General arrangement
    
+                
    Specification
    
+                
+                    
+                        
+                            
+                            
+                            
+                            
+                            
+                        
+                    
+                    
+                        
+                            
+                                                        
+                                                        
+                                                        
+                                                        
+                        
+                        
+                            
+                                                        
+                                                        
+                                                        
+                                                        
+                        
+                        
+                            
+                                                        
+                                                        
+                                                        
+                                                        
+                        
+                        
+                    
+                
+                
    IDLenghtBreadthHeightDeadweight
    Block0 m m m kg
    Block0 m m m kg
    Block2 m m m kg
    
+                
    Calculations
    
+                
+                    
+                        
+                        
+                        
+                    
+                    
+                        
+                             
+                             
+                             
+                        
+                        
+                             
+                             
+                             
+                        
+                        
+                             
+                             
+                             
+                        
+                    
+                
+                
    VariableHand valueCalculated value
    Deadweight (kg)
+                            
    Total area (m2)
+                            
    Total Volume (m3)
+                            
    
+
 			
    
         
     
@@ -342,11 +432,23 @@ 
    Comparation
    
 			controls.update();
 			animate();
             
+            //Ship specification
+            
+            let BOA = ship.structure.hull.attributes.BOA
+            let depth = ship.structure.hull.attributes.Depth
+            let draftdesign = ship.designState.calculationParameters.Draft_design
+            let speeddesign = ship.designState.calculationParameters.speed
+            
+            document.getElementById("LOA").innerHTML = LOA;
+            document.getElementById("BOA").innerHTML = BOA;
+            document.getElementById("depth").innerHTML = depth;
+            document.getElementById("draftdesign").innerHTML = draftdesign;
+            document.getElementById("speeddesign").innerHTML = speeddesign;
+            
             //calculate stability
             let a = ship.calculateStability();
-            let GMc = a.GM.toFixed(3);
-            document.getElementById("GMc").innerHTML = GMc;
-            document.getElementById("GM2").innerHTML = GMc;
+            let GM2 = a.GM.toFixed(3);
+            document.getElementById("GM2").innerHTML = GM2;    
             
             //hand values
             let disp = 100;
@@ -369,8 +471,7 @@ 
    Comparation
    
             let heel = 0;
             let maxmom = 40.072;
             let shear = 13.327;
-            let lightweight = 30000;
-            let deadweight = 70000;
+            let deadweight = 85000;
             
             document.getElementById("disp").innerHTML = disp.toFixed(3);
             document.getElementById("draftms").innerHTML = draftms.toFixed(3);
@@ -392,26 +493,71 @@ 
    Comparation
    
             document.getElementById("heel").innerHTML = heel.toFixed(3);
             document.getElementById("maxmom").innerHTML = maxmom.toFixed(3);
             document.getElementById("shear").innerHTML = shear.toFixed(3);
-            document.getElementById("lightweight").innerHTML = lightweight.toFixed(3);
             document.getElementById("deadweight").innerHTML = deadweight.toFixed(3);
             
+            //baseObjects
+            let lightweight = ship.structure.hull.attributes.lightweight;
+            
+            let block0 = ship.baseObjects[Object.keys(ship.baseObjects)[0]];
+            let block1 = ship.baseObjects[Object.keys(ship.baseObjects)[1]];
+            let block2 = ship.baseObjects[Object.keys(ship.baseObjects)[2]];
+            
+            let bl0 = block0.boxDimensions.length;
+            let bb0 = block0.boxDimensions.breadth;
+            let bh0 = block0.boxDimensions.height;
+            let bwe0 = block0.weightInformation.lightweight;
+            
+            let bl1 = block1.boxDimensions.length;
+            let bb1 = block1.boxDimensions.breadth;
+            let bh1 = block1.boxDimensions.height;
+            let bwe1 = block1.weightInformation.lightweight;
+            
+            let bl2 = block2.boxDimensions.length;
+            let bb2 = block2.boxDimensions.breadth;
+            let bh2 = block2.boxDimensions.height;
+            let bwe2 = block2.weightInformation.lightweight;
+                        
+            let deadweightc = bwe0 + bwe1 + bwe2;
+            
+            let barea = bl0 * bb0 + bl1 * bb1 + bl2 * bb2;
+            let bvol = bl0 * bb0 * bh0 + bl1 * bb1 * bh1 + bl2 * bb2 * bh2;
+                        
+            document.getElementById("bl0").innerHTML = bl0;
+            document.getElementById("bb0").innerHTML = bb0;
+            document.getElementById("bh0").innerHTML = bh0;
+            document.getElementById("bwe0").innerHTML = bwe0;
+            
+            document.getElementById("bl1").innerHTML = bl1;
+            document.getElementById("bb1").innerHTML = bb1;
+            document.getElementById("bh1").innerHTML = bh1;
+            document.getElementById("bwe1").innerHTML = bwe1;
+            
+            document.getElementById("bl2").innerHTML = bl2;
+            document.getElementById("bb2").innerHTML = bb2;
+            document.getElementById("bh2").innerHTML = bh2;
+            document.getElementById("bwe2").innerHTML = bwe2;
+            
+            document.getElementById("barea").innerHTML = barea;
+            document.getElementById("bvol").innerHTML = bvol;
+            
+            //baseObjects calculated values
+            
+            let bareac = 84;
+            let bvolc = 240;
+            document.getElementById("bareac").innerHTML = bareac;
+            document.getElementById("bvolc").innerHTML = bvolc;
+                    
             //calculated values
             let KBc = a.KB.toFixed(3);
             let KGc = a.KG.toFixed(3);
             let BMc = a.BM.toFixed(3);
             let draftc = ship.designState.calculationParameters.Draft_design;
-            
-            //baseObjects not sure why is not working
-            let lightweightc = ship.structure.hull.attributes.lightweight + 15000;
-            let deadweightc = 14000 + 56000;
-            
-            
             let lwlc = ship.structure.hull.calculateAttributesAtDraft(draftc).LWL;
             let bwlc = ship.structure.hull.calculateAttributesAtDraft(draftc).BWL;
             let Awpc = ship.structure.hull.calculateAttributesAtDraft(draftc).Awp;
             let cbc = ship.structure.hull.calculateAttributesAtDraft(draftc).Cb;
             let dispc = lwlc * bwlc * cbc * draftc;
-            
+
             document.getElementById("KBc").innerHTML = KBc;
             document.getElementById("KGc").innerHTML = KGc;
             document.getElementById("BMc").innerHTML = BMc;
@@ -420,7 +566,7 @@ 
    Comparation
    
             document.getElementById("bwlc").innerHTML = bwlc.toFixed(3);
             document.getElementById("Awpc").innerHTML = Awpc.toFixed(3);
             document.getElementById("cbc").innerHTML = cbc.toFixed(3);
-            document.getElementById("lightweightc").innerHTML = lightweightc.toFixed(3);
+            document.getElementById("lightweight").innerHTML = lightweight.toFixed(3);
             document.getElementById("deadweightc").innerHTML = deadweightc.toFixed(3);
             document.getElementById("dispc").innerHTML = dispc.toFixed(3);
             
@@ -428,11 +574,13 @@ 
    Comparation
    
             let KBd = (1-KB/KBc)*100;
             let KGd = (1-KG/KGc)*100;
             let BMd = (1-BM/BMc)*100;
-            let GMd = (1-GM/GMc)*100;
+            let GMd = (1-GM/GM2)*100;
             let draftd = (1-draftms/draftc)*100;
             let lwld = (1-lwl/lwlc)*100;
             let bwld = (1-bwl/bwlc)*100;
             let awpd = (1-Awp/Awpc)*100;
+            let cbd = (1-Cb/cbc)*100;
+            let dispd = (1-disp/dispc)*100;
             
             document.getElementById("KBd").innerHTML = KBd.toFixed(3);
             document.getElementById("KGd").innerHTML = KGd.toFixed(3);
@@ -442,6 +590,8 @@ 
    Comparation
    
             document.getElementById("lwld").innerHTML = lwld.toFixed(3);
             document.getElementById("bwld").innerHTML = bwld.toFixed(3);
             document.getElementById("awpd").innerHTML = awpd.toFixed(3);
+            document.getElementById("cbd").innerHTML = cbd.toFixed(3);
+            document.getElementById("dispd").innerHTML = dispd.toFixed(3);
             
             
             //console.log(ship.structure.hull.calculateAttributesAtDraft(draftc))

From 4602b8ded2836e7aedc1014d3a1bbc6890154fd0 Mon Sep 17 00:00:00 2001
From: EliasHasle 
Date: Fri, 10 Nov 2017 12:24:25 +0100
Subject: [PATCH 20/20] Corrected all hull calculations and stability
 calculations, and added Cp.

---
 build/Vessel.js                        | 186 +++++++++++++------------
 build/makeBuild.py                     |   3 +-
 examples/blockCase_compare_json.html   |   2 +
 examples/js/Ship3D.js                  |   4 +-
 source/CoreClasses/Hull.js             |  14 +-
 source/CoreClasses/Ship.js             |  22 +--
 source/functions/interpolation.js      |  31 +++--
 source/functions/volumeCalculations.js |  98 ++++++-------
 8 files changed, 190 insertions(+), 170 deletions(-)

diff --git a/build/Vessel.js b/build/Vessel.js
index 1f928c5..a6d8494 100644
--- a/build/Vessel.js
+++ b/build/Vessel.js
@@ -1,4 +1,4 @@
-//Vessel.js library, built 2017-11-03 21:58:34.624540, Checksum: 66754104a0f843a0cff421b5aa24f15e
+//Vessel.js library, built 2017-11-10 12:21:12.376049, Checksum: 3162da26a8bf9803d2eabf6c27b72fb4
 /*
 Import like this in HTML:
 
@@ -116,11 +116,13 @@ function linearFromArrays(xx, yy, x) {
 	return lerp(yy[index], yy[index+1], mu);	
 }
 
+//Source: https://en.wikipedia.org/wiki/Bilinear_interpolation
+//(I have used other sources too)
 function bilinearUnitSquareCoeffs(z00, z01, z10, z11) {
-	let a00 = z00;
-	let a10 = z10-z00;
-	let a01 = z01-z00;
-	let a11 = z11+z00-z01-z10;
+	let a00 = z00;				//mux=muy=0
+	let a10 = z10-z00;			//mux=1, muy=0
+	let a01 = z01-z00;			//mux=0, muy=1
+	let a11 = z11+z00-z01-z10;	//mux=muy=1
 	return [a00,a10,a01,a11];
 }
 
@@ -131,7 +133,7 @@ function bilinearUnitSquare(z00, z01, z10, z11, mux, muy) {
 
 //Find coefficients for 1, x, y, xy.
 //This doesn't yet handle zero-lengths well.
-function bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22) {
+function bilinearCoeffs(x1, x2, y1, y2, z00, z01, z10, z11) {
 	let X = (x2-x1);
 	let Y = (y2-y1);
 	
@@ -143,13 +145,13 @@ function bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22) {
 	let Ainv = 1/(X*Y);
 
 	//constant coeff:
-	let b00 = Ainv*(z11*x2*y2 - z21*x1*y2 - z12*x2*y1 + z22*x1*y1);
+	let b00 = Ainv*(z00*x2*y2 - z10*x1*y2 - z01*x2*y1 + z11*x1*y1);
 	//x coeff:
-	let b10 = Ainv*(-z11*y2 + z21*y2 + z12*y1 - z22*y1);
+	let b10 = Ainv*(-z00*y2 + z10*y2 + z01*y1 - z11*y1);
 	//y coeff:
-	let b01 = Ainv*(-z11*x2 + z21*x1 + z12*x2 -z22*x1);
+	let b01 = Ainv*(-z00*x2 + z10*x1 + z01*x2 -z11*x1);
 	//xy coeff:
-	let b11 = Ainv*(z11-z21-z12+z22);
+	let b11 = Ainv*(z00-z10-z01+z11);
 	
 	return [b00,b10,b01,b11];
 }
@@ -160,11 +162,16 @@ function bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22) {
 function bilinear(x1, x2, y1, y2, z11, z12, z21, z22, x, y) {
 	let [b00, b10, b01, b11] = 
 		bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22);
-	return b00 + b10*x + b01*y + b11*x*y;
-	//The following is supposed to be equivalent. Maybe I should compare, to make sure that the current calculation is correct.
+	let fromCoeffs = b00 + b10*x + b01*y + b11*x*y;
+
+	//The following is supposed to be equivalent. Some tests yielding identical results (and no tests so far yielding different results) suggest that the calculations are in fact equivalent.
 	/*let mux = (x-x1)/(x2-x1);
 	let muy = (y-y1)/(y2-y1);
-	return bilinearUnitSquare(z11, z12, z21, z22, mux, muy);*/
+	let fromUnitSquare = bilinearUnitSquare(z11, z12, z21, z22, mux, muy);
+	
+	console.log("fromCoeffs=", fromCoeffs, ", fromUnitSquare=", fromUnitSquare);*/
+	
+	return fromCoeffs;
 }
 //@EliasHasle
 
@@ -354,67 +361,67 @@ function signedTriangleArea(v1,v2,v3) {
 	return A;
 }//@EliasHasle
 
-//This one is broken, apparently.
-//It returns negative xc and yc when the z values are negative, 
-//even though all x and y values are positive.
-//I am currently doing some tests here of an (over-)simplified calculation
-//The results so far indicate that, for the prism hull, the results are almost identical, except that with the simple calculation the center of volume is almost right (but wrong enough to disqualify such a simple calculation). The bug causing very wrong (too big) submerged volume is likely to be found elsewhere.
-/*Note that the coordinate system used here has xy as a grid, with z as heights on the grid, but in the intended application, which is calculations on transverse hull offsets, z corresponds to the vessel y axis, and y corresponds to the vessel z axis. In the application, the conversion between coordinate systems must be taken care of appropriately.*/
+//I have been doing some tests here of a simplified calculation.
+//The results so far indicate that, for the prism hull, the results are almost identical, except that with the simple calculation the center of volume is almost right (but wrong enough to disqualify such a simple calculation).
+/*Note that the coordinate system used here has xy as a grid, with z as heights on the grid, but in the intended application, which is calculations on transverse hull offsets, this z corresponds to the vessel y axis, and y corresponds to the vessel z axis. In any application of this function, the conversion between coordinate systems must be taken care of appropriately.*/
 											  // xy
-function patchColumnCalculation(x1, x2, y1, y2, z11, z12, z21, z22) {
-	//DEBUG START:
-	//Simpler approximate calculation of volume:
-	let z = 0.25*(z11+z12+z21+z22); //"average height"
-	let V = Math.abs((x2-x1)*(y2-y1)*z); //base area times "average height"
-
-	//Very approximate center of volume
-	//(does not account for different signs on z values,
-	//but that should be OK for hull offsets)
-	let xc = (x1*(z11+z12)+x2*(z21+z22))/((z11+z12+z21+z22) || 1);
-	let yc = (y1*(z11+z21)+y2*(z12+z22))/((z11+z12+z21+z22) || 1);
-	let zc = 0.5*z;
-	
-	//Simple triangle average approximation for area
-	let As = elementArea(
-		{x: x1, y: y1, z: z11},
-		{x: x1, y: y2, z: z12},
-		{x: x2, y: y1, z: z21},
-		{x: x2, y: y2, z: z22});
-	
-	return {As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
-	
-	//DEBUG END
-	
-	//Calculation based on a bilinear patch:
-
-	// let X = x2-x1;
-	// let Y = y2-y1;
-	// let [a00, a10, a01, a11] = bilinearUnitSquareCoeffs(z11, z12, z21, z22);
+function patchColumnCalculation(x1, x2, y1, y2, z00, z01, z10, z11) {
+	//VOLUME:
+	//Analysis based on a bilinear patch:
 	// /*
 	// From here I call mux for x, and muy for y.
 	// Integral over unit square:
 	// INT[x from 0 to 1, INT[y from 0 to 1, (a00 + a10*x + a01*y + a11*x*y) dy] dx]
 	// = INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x) dx]
 	// = a00 + 0.5*a10 + 0.5*a01 + 0.25*a11
+	// Note that by expanding a00,a10,a01,a11, it is demonstrated that this (rather unsurprisingly) is exactly equivalent to taking the average z offset of the control points.
 	// */
-	// let Ab = X*Y; //area of base of patch column
-	// let zAvg = (a00 + 0.5*a10 + 0.5*a01 + 0.25*a11);
-	// let V = Math.abs(Ab*zAvg); //new: absolute value
-	// let zc = 0.5*zAvg;
+	 let X = x2-x1;
+	 let Y = y2-y1;
+	 let Ab = X*Y; //area of base of patch column
+	 //let zAvg = (a00 + 0.5*a10 + 0.5*a01 + 0.25*a11);
+	 let zAvg = 0.25*(z00+z01+z10+z11); //equivalent
+	 let V = Math.abs(Ab*zAvg); //works
+
+	//CENTER OF VOLUME
+	 let zc = 0.5*zAvg;
+	 
+	//Very approximate center of volume
+	//(does not account for different signs on z values,
+	//but that should be OK for hull offsets)
+	//let xc = (x1*(z00+z01)+x2*(z10+z11))/((z00+z01+z10+z11) || 1);
+	//let yc = (y1*(z00+z10)+y2*(z01+z11))/((z00+z01+z10+z11) || 1);
+	
 	// /*
-	// To find xc, I need to integrate x*z over the unit square, and scale and translate to ship coordinates afterwards:
-	// INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x)*x dx]
-	// = 0.5*a00 + a10/3 + 0.25*a01 + a11/6
-	// Scale and translate:*/
-	// let xc = x1 + X*(0.5*a00 + a10/3 + 0.25*a01 + a11/6);
+	// To find xc properly, I need to integrate x*z over the unit square, divide by zAvg(?) and scale and translate to ship coordinates afterwards:
+	// INT[x from 0 to 1, INT[y from 0 to 1, x*(a00 + a10*x + a01*y + a11*x*y) dy] dx] =
+	// INT[x from 0 to 1, INT[y from 0 to 1, (a00*x + a10*x^2 + a01*xy + a11*x^2*y) dy] dx] =
+	// INT[x from 0 to 1, (a00*x + a10*x^2 + 0.5*a01*x + 0.5*a11*x^2) dx]
+	// = (0.5*a00 + a10/3 + 0.25*a01 + a11/6)
+	//Trying to expand the coeffs to original z offsets:
+	// = (0.5*z00 + (z10-z00)/3 + 0.25*(z01-z00) + (z00+z00-z01-z10)/6)
+	// = ((1/12)*z00 + (1/6)*z10 + (1/12)*z01 + (1/6)*z00)
+	//Divide by zAvg to get muxc, then scale and translate to xc.
+	let xc = x1+X*(((1/12)*z00 + (1/6)*z10 + (1/12)*z01 + (1/6)*z11) / (zAvg || 1));
+	//console.log("x1=%.2f, X=%.2f, muxc = %.2f", x1, X, (((1/12)*z00 + (1/6)*z10 + (1/12)*z01 + (1/6)*z11) / (zAvg || 1)));
+	//Similar for yc (modified symmetrically)
+	let yc = y1+Y*(((1/12)*z00 + (1/12)*z10 + (1/6)*z01 + (1/6)*z11) / (zAvg || 1));
+	let [a00, a10, a01, a11] = bilinearUnitSquareCoeffs(z00, z01, z10, z11);
 	
-	// //Similar for yc:
-	// let yc = y1 + Y*(0.5*a00 + 0.25*a10 + a01/3 + a11/6);
+	//console.log("Patch column Cv = (%.2f, %.2f, %.2f)", xc,yc,zc);
 	
-	//new: absolute value (OK?)
-	//let As = Math.abs(bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22));
+	//AREA
+	//These two methods give very similar results, within about 1% difference for the fishing boat hull (used in PX121.json).
+	//Simple triangle average approximation for area (works)
+	/*let As = elementArea(
+		{x: x1, y: y1, z: z00},
+		{x: x1, y: y2, z: z01},
+		{x: x2, y: y1, z: z10},
+		{x: x2, y: y2, z: z11});*/
+	//Bilinear area calculation. Works too, but is currently numerical, and quite complex (which means it is bug-prone and hard to maintain). But it is more exact, even with just a few segments for numerical integration (the last, optional, parameter)
+	let As = Math.abs(bilinearArea(x1, x2, y1, y2, z00, z01, z10, z11, 10));
 	
-	// return {Ab: Ab, As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
+	return {Ab: Ab, As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
 }
 
 //Input: array of objects with calculation results for elements.
@@ -709,6 +716,8 @@ Object.assign(Ship.prototype, {
 
 		return specification;
 	},
+	//This should probably be separated in lightweight and deadweight
+	//Then this function should be replaced by a getDisplacement
 	getWeight: function(shipState) {
 		shipState = shipState || this.designState;
 
@@ -729,6 +738,7 @@ Object.assign(Ship.prototype, {
 		console.info("Calculated weight object: ", W);
 		return W;
 	},
+	//This should just take displacement as parameter instead. (later, soon)
 	calculateDraft(shipState, epsilon=0.001) {
 		let w = this.getWeight(shipState);
 		let M = w.mass;
@@ -752,18 +762,15 @@ Object.assign(Ship.prototype, {
         let T = this.calculateDraft(shipState);
         let ha = this.structure.hull.calculateAttributesAtDraft(T);
         let vol = ha.Vs;
-        if (vol === 0){
-            let Lwl = this.designState.calculationParameters.LWL_design;
-            let B = this.structure.hull.attributes.BOA;
-            let cb = this.designState.calculationParameters.Cb_design;
-            vol = Lwl * B * T * cb;
-        }
         let KG = this.getWeight(shipState).cg.z;
-        let I = ha.Iywp;
-        let KB = 0.52 * T;
-        let BM = I / vol;
-        let GM = KB + BM - KG;
-        return {GM, KB, BM, KG};
+		let Ix = ha.Ixwp;
+        let Iy = ha.Iywp;
+        let KB = ha.Cv.z;
+		let BMT = Ix / vol;
+        let BML = Iy / vol;
+        let GMT = KB + BMT - KG;
+        let GML = KB + BML - KG;
+        return {GMT, GML, GM: T, KB, BMT, BML, BM: BMT, KG};
 	}
 });//@EliasHasle
 
@@ -830,7 +837,7 @@ Object.assign(Structure.prototype, {
 			let zc = d.zFloor+0.5*d.thickness;
 			let yc = d.yCentre;
 			let b = d.breadth;
-			let wlc = this.hull.waterlineCalculation(zc, {minX: d.xAft, maxX: d.xFwd, minY: yc-0.5*b, maxY: yc+0.5*b}, 3);
+			let wlc = this.hull.waterlineCalculation(zc, {minX: d.xAft, maxX: d.xFwd, minY: yc-0.5*b, maxY: yc+0.5*b});
 			components.push({
 				//approximation
 				mass: wlc.Awp*d.thickness*d.density,
@@ -915,8 +922,6 @@ Object.assign(Hull.prototype, {
 	/*
 	Testing new version without nanCorrectionMode parameter, that defaults to setting lower NaNs to 0 and extrapolating highest data entry for upper NaNs (if existant, else set to 0). Inner NaNs will also be set to zero.
 	
-	This does not exactly work perfectly yet. getwaterline(0,3) gives interpolated values, even though the keel level is defined.
-	
 	Input:
 	z: level from bottom of ship (absolute value in meters)
 	
@@ -1043,13 +1048,13 @@ Object.assign(Hull.prototype, {
 	//THIS is a candidate for causing wrong Ix, Iy values.
 	//Much logic that can go wrong.
 									//typically deck bounds
-	waterlineCalculation: function(z, bounds, nanCorrectionMode=3) {
+	waterlineCalculation: function(z, bounds) {
 		let {minX, maxX, minY, maxY} = bounds || {};
 
 		console.group/*Collapsed*/("waterlineCalculation.");
-		console.info("Arguments: z=", z, " Boundaries: ", arguments[1], " NaN correction mode: ", nanCorrectionMode);
+		console.info("Arguments: z=", z, " Boundaries: ", arguments[1]);
 		
-		let wl = this.getWaterline(z, nanCorrectionMode);
+		let wl = this.getWaterline(z);
 		console.info("wl: ", wl); //DEBUG
 
 		let LOA = this.attributes.LOA;
@@ -1177,7 +1182,7 @@ Object.assign(Hull.prototype, {
 		};
 	},
 
-	//NOT DONE YET. Calculates too big Ap, Vs, Cb, and too small As for the test case. For the test, the Ap is exactly BWL*1m larger than it should be (why?). Too high Cb is clearly caused by too big Vs, and big Ap and Vs may have a common cause. The bilinear volume and area calculations have been temporarily replaced with simpler calculations, but this does not seem to help. I expect to find the bug(s) elsewhere.
+	//NOT DONE YET. Many bugs are fixed, but the volume center calculation is broken. The bilinear volume and area calculations have been temporarily replaced with simpler (and worse) calculations, and at least the volume and volume center calculations should be revived soon.
 	//Important: calculateAttributesAtDraft takes one mandatory parameter T. (The function defined here is immediately called during construction of the prototype, and returns the proper function.)
 	calculateAttributesAtDraft: function() {
 		function levelCalculation(hull,
@@ -1197,7 +1202,7 @@ Object.assign(Hull.prototype, {
 				Cv: {x:0, y:0, z:0}
 			}) {
 			
-			let wlc = hull.waterlineCalculation(z,{},3);
+			let wlc = hull.waterlineCalculation(z,{});
 			let lev = {};
 			Object.assign(lev, wlc);
 			//Projected area calculation (approximate):
@@ -1242,8 +1247,8 @@ Object.assign(Hull.prototype, {
 			//Many possibilities for getting the coordinate systems wrong.
 			let calculations = [];
 			let sts = hull.halfBreadths.stations.map(st=>st*hull.attributes.LOA);
-			let wl = hull.getWaterline(z,3);
-			let prwl = hull.getWaterline(prev.z,3);
+			let wl = hull.getWaterline(z);
+			let prwl = hull.getWaterline(prev.z);
 			for (let j = 0; j < sts.length-1; j++) {
 				let port = 
 					patchColumnCalculation(sts[j], sts[j+1], prev.z, z, -prwl[j], -wl[j], -prwl[j+1], -wl[j+1]);
@@ -1274,6 +1279,7 @@ Object.assign(Hull.prototype, {
 					), 1/(lev.Vs || 2));
 			
 			lev.Cb = lev.Vs/lev.Vbb;
+			lev.Cp = lev.Vs/(lev.Ap*(lev.maxX-lev.minX));
 			
 			return lev;
 		}
@@ -1311,7 +1317,7 @@ Object.assign(Hull.prototype, {
 			
 			//Filter and rename for output
 			return {
-				xcwp: lc.xc,
+				xcwp: lc.xc, //water plane values
 				ycwp: lc.yc,
 				Awp: lc.Awp,
 				Ixwp: lc.Ix,
@@ -1324,12 +1330,13 @@ Object.assign(Hull.prototype, {
 				LWL: lc.LWL,
 				LBP: lc.LBP,
 				BWL: lc.BWL,
-				Ap: lc.Ap,
+				Ap: lc.Ap, //projected area in length direction
+				Cp: lc.Cp, //prismatic coefficient
 				//Vbb: lc.Vbb,
-				Vs: lc.Vs,
+				Vs: lc.Vs, //volume of submerged part of the hull
 				Cb: lc.Cb,
-				As: lc.As,
-				Cv: lc.Cv			
+				As: lc.As, //wetted area
+				Cv: lc.Cv //center of buoyancy
 			}
 		};
 	}()
@@ -1749,7 +1756,8 @@ Object.assign(Vessel, {
 	loadShip: loadShip,
 	downloadShip: downloadShip,
         f: {
-            linearFromArrays: linearFromArrays
+            linearFromArrays: linearFromArrays,
+            bilinear: bilinear
         }
 });
 })();
diff --git a/build/makeBuild.py b/build/makeBuild.py
index 86b0cfb..67b67a5 100644
--- a/build/makeBuild.py
+++ b/build/makeBuild.py
@@ -43,7 +43,8 @@
 	loadShip: loadShip,
 	downloadShip: downloadShip,
         f: {
-            linearFromArrays: linearFromArrays
+            linearFromArrays: linearFromArrays,
+            bilinear: bilinear
         }
 });
 })();
diff --git a/examples/blockCase_compare_json.html b/examples/blockCase_compare_json.html
index 8021632..ec75779 100644
--- a/examples/blockCase_compare_json.html
+++ b/examples/blockCase_compare_json.html
@@ -100,6 +100,7 @@
 	let Vs = 0.5*LWL*BWL*T;
 	//let Vbb = BWL*LWL*T;// = LOA*(BOA/Depth)*T^2
 	let Cb = 0.5;// = Vs/Vbb
+	let Cp = 1; // because the hull is a prism
 	let Ap = 0.5*BWL*T;
 	let As = 2*Ap + 2*LWL*Math.sqrt(T**2 + (0.5*BWL)**2);
 	let Awp = LWL*BWL;
@@ -127,6 +128,7 @@
 		Awp: Awp,
 		Cwp: Cwp,
 		Cb: Cb,
+		Cp: Cp,
 		Cv: Cv,/*{
 			x: LCBx,
 			y: undefined,
diff --git a/examples/js/Ship3D.js b/examples/js/Ship3D.js
index 67cd636..1be5b4f 100644
--- a/examples/js/Ship3D.js
+++ b/examples/js/Ship3D.js
@@ -154,8 +154,8 @@ function Ship3D(ship, stlPath) {
 		console.log("d.zFloor=%.1f", d.zFloor); //DEBUG
 		let zHigh = d.zFloor;
 		let zLow = d.zFloor-d.thickness;
-		let wlHigh = ship.structure.hull.getWaterline(zHigh, 2);
-		let wlLow = ship.structure.hull.getWaterline(zLow, 2);
+		let wlHigh = ship.structure.hull.getWaterline(zHigh);
+		let wlLow = ship.structure.hull.getWaterline(zLow);
 		let pos = deckGeom.getAttribute("position");
 		let pa = pos.array;
 		for (let j = 0; j < stss.length+1; j++) {
diff --git a/source/CoreClasses/Hull.js b/source/CoreClasses/Hull.js
index 7505223..2f9983d 100644
--- a/source/CoreClasses/Hull.js
+++ b/source/CoreClasses/Hull.js
@@ -311,7 +311,7 @@ Object.assign(Hull.prototype, {
 		};
 	},
 
-	//NOT DONE YET. Calculates too big Ap, Vs, Cb, and too small As for the test case. For the test, the Ap is exactly BWL*1m larger than it should be (why?). Too high Cb is clearly caused by too big Vs, and big Ap and Vs may have a common cause. The bilinear volume and area calculations have been temporarily replaced with simpler calculations, but this does not seem to help. I expect to find the bug(s) elsewhere.
+	//NOT DONE YET. Many bugs are fixed, but the volume center calculation is broken. The bilinear volume and area calculations have been temporarily replaced with simpler (and worse) calculations, and at least the volume and volume center calculations should be revived soon.
 	//Important: calculateAttributesAtDraft takes one mandatory parameter T. (The function defined here is immediately called during construction of the prototype, and returns the proper function.)
 	calculateAttributesAtDraft: function() {
 		function levelCalculation(hull,
@@ -408,6 +408,7 @@ Object.assign(Hull.prototype, {
 					), 1/(lev.Vs || 2));
 			
 			lev.Cb = lev.Vs/lev.Vbb;
+			lev.Cp = lev.Vs/(lev.Ap*(lev.maxX-lev.minX));
 			
 			return lev;
 		}
@@ -445,7 +446,7 @@ Object.assign(Hull.prototype, {
 			
 			//Filter and rename for output
 			return {
-				xcwp: lc.xc,
+				xcwp: lc.xc, //water plane values
 				ycwp: lc.yc,
 				Awp: lc.Awp,
 				Ixwp: lc.Ix,
@@ -458,12 +459,13 @@ Object.assign(Hull.prototype, {
 				LWL: lc.LWL,
 				LBP: lc.LBP,
 				BWL: lc.BWL,
-				Ap: lc.Ap,
+				Ap: lc.Ap, //projected area in length direction
+				Cp: lc.Cp, //prismatic coefficient
 				//Vbb: lc.Vbb,
-				Vs: lc.Vs,
+				Vs: lc.Vs, //volume of submerged part of the hull
 				Cb: lc.Cb,
-				As: lc.As,
-				Cv: lc.Cv			
+				As: lc.As, //wetted area
+				Cv: lc.Cv //center of buoyancy
 			}
 		};
 	}()
diff --git a/source/CoreClasses/Ship.js b/source/CoreClasses/Ship.js
index b8f1131..f83854f 100644
--- a/source/CoreClasses/Ship.js
+++ b/source/CoreClasses/Ship.js
@@ -45,6 +45,8 @@ Object.assign(Ship.prototype, {
 
 		return specification;
 	},
+	//This should probably be separated in lightweight and deadweight
+	//Then this function should be replaced by a getDisplacement
 	getWeight: function(shipState) {
 		shipState = shipState || this.designState;
 
@@ -65,6 +67,7 @@ Object.assign(Ship.prototype, {
 		console.info("Calculated weight object: ", W);
 		return W;
 	},
+	//This should just take displacement as parameter instead. (later, soon)
 	calculateDraft(shipState, epsilon=0.001) {
 		let w = this.getWeight(shipState);
 		let M = w.mass;
@@ -88,17 +91,14 @@ Object.assign(Ship.prototype, {
         let T = this.calculateDraft(shipState);
         let ha = this.structure.hull.calculateAttributesAtDraft(T);
         let vol = ha.Vs;
-        if (vol === 0){
-            let Lwl = this.designState.calculationParameters.LWL_design;
-            let B = this.structure.hull.attributes.BOA;
-            let cb = this.designState.calculationParameters.Cb_design;
-            vol = Lwl * B * T * cb;
-        }
         let KG = this.getWeight(shipState).cg.z;
-        let I = ha.Iywp;
-        let KB = 0.52 * T;
-        let BM = I / vol;
-        let GM = KB + BM - KG;
-        return {GM, KB, BM, KG};
+		let Ix = ha.Ixwp;
+        let Iy = ha.Iywp;
+        let KB = ha.Cv.z;
+		let BMT = Ix / vol;
+        let BML = Iy / vol;
+        let GMT = KB + BMT - KG;
+        let GML = KB + BML - KG;
+        return {GMT, GML, GM: T, KB, BMT, BML, BM: BMT, KG};
 	}
 });
\ No newline at end of file
diff --git a/source/functions/interpolation.js b/source/functions/interpolation.js
index 3101738..abefab7 100644
--- a/source/functions/interpolation.js
+++ b/source/functions/interpolation.js
@@ -44,11 +44,13 @@ function linearFromArrays(xx, yy, x) {
 	return lerp(yy[index], yy[index+1], mu);	
 }
 
+//Source: https://en.wikipedia.org/wiki/Bilinear_interpolation
+//(I have used other sources too)
 function bilinearUnitSquareCoeffs(z00, z01, z10, z11) {
-	let a00 = z00;
-	let a10 = z10-z00;
-	let a01 = z01-z00;
-	let a11 = z11+z00-z01-z10;
+	let a00 = z00;				//mux=muy=0
+	let a10 = z10-z00;			//mux=1, muy=0
+	let a01 = z01-z00;			//mux=0, muy=1
+	let a11 = z11+z00-z01-z10;	//mux=muy=1
 	return [a00,a10,a01,a11];
 }
 
@@ -59,7 +61,7 @@ function bilinearUnitSquare(z00, z01, z10, z11, mux, muy) {
 
 //Find coefficients for 1, x, y, xy.
 //This doesn't yet handle zero-lengths well.
-function bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22) {
+function bilinearCoeffs(x1, x2, y1, y2, z00, z01, z10, z11) {
 	let X = (x2-x1);
 	let Y = (y2-y1);
 	
@@ -71,13 +73,13 @@ function bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22) {
 	let Ainv = 1/(X*Y);
 
 	//constant coeff:
-	let b00 = Ainv*(z11*x2*y2 - z21*x1*y2 - z12*x2*y1 + z22*x1*y1);
+	let b00 = Ainv*(z00*x2*y2 - z10*x1*y2 - z01*x2*y1 + z11*x1*y1);
 	//x coeff:
-	let b10 = Ainv*(-z11*y2 + z21*y2 + z12*y1 - z22*y1);
+	let b10 = Ainv*(-z00*y2 + z10*y2 + z01*y1 - z11*y1);
 	//y coeff:
-	let b01 = Ainv*(-z11*x2 + z21*x1 + z12*x2 -z22*x1);
+	let b01 = Ainv*(-z00*x2 + z10*x1 + z01*x2 -z11*x1);
 	//xy coeff:
-	let b11 = Ainv*(z11-z21-z12+z22);
+	let b11 = Ainv*(z00-z10-z01+z11);
 	
 	return [b00,b10,b01,b11];
 }
@@ -88,9 +90,14 @@ function bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22) {
 function bilinear(x1, x2, y1, y2, z11, z12, z21, z22, x, y) {
 	let [b00, b10, b01, b11] = 
 		bilinearCoeffs(x1, x2, y1, y2, z11, z12, z21, z22);
-	return b00 + b10*x + b01*y + b11*x*y;
-	//The following is supposed to be equivalent. Maybe I should compare, to make sure that the current calculation is correct.
+	let fromCoeffs = b00 + b10*x + b01*y + b11*x*y;
+
+	//The following is supposed to be equivalent. Some tests yielding identical results (and no tests so far yielding different results) suggest that the calculations are in fact equivalent.
 	/*let mux = (x-x1)/(x2-x1);
 	let muy = (y-y1)/(y2-y1);
-	return bilinearUnitSquare(z11, z12, z21, z22, mux, muy);*/
+	let fromUnitSquare = bilinearUnitSquare(z11, z12, z21, z22, mux, muy);
+	
+	console.log("fromCoeffs=", fromCoeffs, ", fromUnitSquare=", fromUnitSquare);*/
+	
+	return fromCoeffs;
 }
diff --git a/source/functions/volumeCalculations.js b/source/functions/volumeCalculations.js
index 1f8b321..e8891e5 100644
--- a/source/functions/volumeCalculations.js
+++ b/source/functions/volumeCalculations.js
@@ -1,66 +1,66 @@
 //@EliasHasle
 
-//This one is broken, apparently.
-//It returns negative xc and yc when the z values are negative, 
-//even though all x and y values are positive.
-//I am currently doing some tests here of an (over-)simplified calculation
-//The results so far indicate that, for the prism hull, the results are almost identical, except that with the simple calculation the center of volume is almost right (but wrong enough to disqualify such a simple calculation). The bug causing very wrong (too big) submerged volume is likely to be found elsewhere.
-/*Note that the coordinate system used here has xy as a grid, with z as heights on the grid, but in the intended application, which is calculations on transverse hull offsets, z corresponds to the vessel y axis, and y corresponds to the vessel z axis. In the application, the conversion between coordinate systems must be taken care of appropriately.*/
+//I have been doing some tests here of a simplified calculation.
+//The results so far indicate that, for the prism hull, the results are almost identical, except that with the simple calculation the center of volume is almost right (but wrong enough to disqualify such a simple calculation).
+/*Note that the coordinate system used here has xy as a grid, with z as heights on the grid, but in the intended application, which is calculations on transverse hull offsets, this z corresponds to the vessel y axis, and y corresponds to the vessel z axis. In any application of this function, the conversion between coordinate systems must be taken care of appropriately.*/
 											  // xy
-function patchColumnCalculation(x1, x2, y1, y2, z11, z12, z21, z22) {
-	//DEBUG START:
-	//Simpler approximate calculation of volume:
-	let z = 0.25*(z11+z12+z21+z22); //"average height"
-	let V = Math.abs((x2-x1)*(y2-y1)*z); //base area times "average height"
-
-	//Very approximate center of volume
-	//(does not account for different signs on z values,
-	//but that should be OK for hull offsets)
-	let xc = (x1*(z11+z12)+x2*(z21+z22))/((z11+z12+z21+z22) || 1);
-	let yc = (y1*(z11+z21)+y2*(z12+z22))/((z11+z12+z21+z22) || 1);
-	let zc = 0.5*z;
-	
-	//Simple triangle average approximation for area
-	let As = elementArea(
-		{x: x1, y: y1, z: z11},
-		{x: x1, y: y2, z: z12},
-		{x: x2, y: y1, z: z21},
-		{x: x2, y: y2, z: z22});
-	
-	return {As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
-	
-	//DEBUG END
-	
-	//Calculation based on a bilinear patch:
-
-	// let X = x2-x1;
-	// let Y = y2-y1;
-	// let [a00, a10, a01, a11] = bilinearUnitSquareCoeffs(z11, z12, z21, z22);
+function patchColumnCalculation(x1, x2, y1, y2, z00, z01, z10, z11) {
+	//VOLUME:
+	//Analysis based on a bilinear patch:
 	// /*
 	// From here I call mux for x, and muy for y.
 	// Integral over unit square:
 	// INT[x from 0 to 1, INT[y from 0 to 1, (a00 + a10*x + a01*y + a11*x*y) dy] dx]
 	// = INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x) dx]
 	// = a00 + 0.5*a10 + 0.5*a01 + 0.25*a11
+	// Note that by expanding a00,a10,a01,a11, it is demonstrated that this (rather unsurprisingly) is exactly equivalent to taking the average z offset of the control points.
 	// */
-	// let Ab = X*Y; //area of base of patch column
-	// let zAvg = (a00 + 0.5*a10 + 0.5*a01 + 0.25*a11);
-	// let V = Math.abs(Ab*zAvg); //new: absolute value
-	// let zc = 0.5*zAvg;
+	 let X = x2-x1;
+	 let Y = y2-y1;
+	 let Ab = X*Y; //area of base of patch column
+	 //let zAvg = (a00 + 0.5*a10 + 0.5*a01 + 0.25*a11);
+	 let zAvg = 0.25*(z00+z01+z10+z11); //equivalent
+	 let V = Math.abs(Ab*zAvg); //works
+
+	//CENTER OF VOLUME
+	 let zc = 0.5*zAvg;
+	 
+	//Very approximate center of volume
+	//(does not account for different signs on z values,
+	//but that should be OK for hull offsets)
+	//let xc = (x1*(z00+z01)+x2*(z10+z11))/((z00+z01+z10+z11) || 1);
+	//let yc = (y1*(z00+z10)+y2*(z01+z11))/((z00+z01+z10+z11) || 1);
+	
 	// /*
-	// To find xc, I need to integrate x*z over the unit square, and scale and translate to ship coordinates afterwards:
-	// INT[x from 0 to 1, (a00+a10*x+0.5*a01+0.5*a11*x)*x dx]
-	// = 0.5*a00 + a10/3 + 0.25*a01 + a11/6
-	// Scale and translate:*/
-	// let xc = x1 + X*(0.5*a00 + a10/3 + 0.25*a01 + a11/6);
+	// To find xc properly, I need to integrate x*z over the unit square, divide by zAvg(?) and scale and translate to ship coordinates afterwards:
+	// INT[x from 0 to 1, INT[y from 0 to 1, x*(a00 + a10*x + a01*y + a11*x*y) dy] dx] =
+	// INT[x from 0 to 1, INT[y from 0 to 1, (a00*x + a10*x^2 + a01*xy + a11*x^2*y) dy] dx] =
+	// INT[x from 0 to 1, (a00*x + a10*x^2 + 0.5*a01*x + 0.5*a11*x^2) dx]
+	// = (0.5*a00 + a10/3 + 0.25*a01 + a11/6)
+	//Trying to expand the coeffs to original z offsets:
+	// = (0.5*z00 + (z10-z00)/3 + 0.25*(z01-z00) + (z00+z00-z01-z10)/6)
+	// = ((1/12)*z00 + (1/6)*z10 + (1/12)*z01 + (1/6)*z00)
+	//Divide by zAvg to get muxc, then scale and translate to xc.
+	let xc = x1+X*(((1/12)*z00 + (1/6)*z10 + (1/12)*z01 + (1/6)*z11) / (zAvg || 1));
+	//console.log("x1=%.2f, X=%.2f, muxc = %.2f", x1, X, (((1/12)*z00 + (1/6)*z10 + (1/12)*z01 + (1/6)*z11) / (zAvg || 1)));
+	//Similar for yc (modified symmetrically)
+	let yc = y1+Y*(((1/12)*z00 + (1/12)*z10 + (1/6)*z01 + (1/6)*z11) / (zAvg || 1));
+	let [a00, a10, a01, a11] = bilinearUnitSquareCoeffs(z00, z01, z10, z11);
 	
-	// //Similar for yc:
-	// let yc = y1 + Y*(0.5*a00 + 0.25*a10 + a01/3 + a11/6);
+	//console.log("Patch column Cv = (%.2f, %.2f, %.2f)", xc,yc,zc);
 	
-	//new: absolute value (OK?)
-	//let As = Math.abs(bilinearArea(x1, x2, y1, y2, z11, z12, z21, z22));
+	//AREA
+	//These two methods give very similar results, within about 1% difference for the fishing boat hull (used in PX121.json).
+	//Simple triangle average approximation for area (works)
+	/*let As = elementArea(
+		{x: x1, y: y1, z: z00},
+		{x: x1, y: y2, z: z01},
+		{x: x2, y: y1, z: z10},
+		{x: x2, y: y2, z: z11});*/
+	//Bilinear area calculation. Works too, but is currently numerical, and quite complex (which means it is bug-prone and hard to maintain). But it is more exact, even with just a few segments for numerical integration (the last, optional, parameter)
+	let As = Math.abs(bilinearArea(x1, x2, y1, y2, z00, z01, z10, z11, 10));
 	
-	// return {Ab: Ab, As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
+	return {Ab: Ab, As: As, V: V, Cv: {x: xc, y: yc, z: zc}};
 }
 
 //Input: array of objects with calculation results for elements.