// zforce.cmd

// Some ways to calculate vertical force.
// zforce4 is the preferred way.

dz := 1e-5 * height // perturbation size

// Single difference, moving just top pad
zforce1 := { old_energy := total_energy;
             height := height + dz;
             recalc;
             new_energy := total_energy - body[1].pressure *
                      (body[1].volume - body[1].target);
             zforce := -(new_energy - old_energy)/dz;
             printf "zforce: %18.15f  (single difference, top pad only)\n",
		zforce;
             height := height - dz; // restore original
             recalc;
           }

// Central difference, moving just top pad
zforce2 := { 
             height := height + dz;
             recalc;
             new_energy_1 := total_energy - body[1].pressure *
                      (body[1].volume - body[1].target);
             height := height - 2*dz;
             recalc;
             new_energy_2 := total_energy - body[1].pressure *
                      (body[1].volume - body[1].target);
             zforce := -(new_energy_1 - new_energy_2)/2/dz;
             printf "zforce: %18.15f  (central difference, top pad only)\n",
	       zforce;
             height := height + dz; // restore original
             recalc;
           }

// Single difference, moving linearly
zforce3 := { old_energy := total_energy;
	     set vertex z z + dz*z/height; // do this before changing height
             height := height + dz;
             recalc;
             new_energy := total_energy - body[1].pressure *
                      (body[1].volume - body[1].target);
             zforce := -(new_energy - old_energy)/dz;
             printf "zforce: %18.15f  (single difference, linear move)\n",
		zforce;
	     set vertex z z - dz*z/height; 
             height := height - dz; // restore original
             recalc;
           }

// Central difference, moving linearly
zforce4 := { 
	     set vertex z z + dz*z/height; // do this before changing height
             height := height + dz;
             recalc;
             new_energy_1 := total_energy - body[1].pressure *
                      (body[1].volume - body[1].target);
	     set vertex z z - 2*dz*z/height; 
             height := height - 2*dz;
             recalc;
             new_energy_2 := total_energy - body[1].pressure *
                      (body[1].volume - body[1].target);
             zforce := -(new_energy_1 - new_energy_2)/2/dz;
             printf "zforce: %18.15f  (central difference, linear move)\n",
	       zforce;
	     set vertex z z + dz*z/height; 
             height := height + dz; // restore original
             recalc;
           }


// Central difference, moving linearly, with re-convergence,
// so don't need volume corrections.
zforce5 := { 
	     set vertex z z + dz*z/height; // do this before changing height
             height := height + dz;
             recalc; hessian; hessian;
             new_energy_1 := total_energy; 
	     set vertex z z - 2*dz*z/height; 
             height := height - 2*dz;
             recalc; hessian; hessian;
             new_energy_2 := total_energy;
             zforce := -(new_energy_1 - new_energy_2)/2/dz;
             printf "zforce: %18.15f  (central difference, linear move, reconverge)\n",
	       zforce;
	     set vertex z z + dz*z/height; 
             height := height + dz; // restore original
             recalc; hessian; hessian;  // reconverge
           }

// Geometrical calculation of forces
// (not accurate since it leaves out effect of facets touching
// top pad at only one vertex)

zforce6 := { // first, pressure on top, accounting for gravity 
	     zforce := sum(facet where on_constraint 2,area) *
               (body[1].pressure - gravity_constant*height*SOLDER_DENSITY);
	     // add tension due to facets pulling on top
	     foreach edge ee where on_constraint 2 do
	     { zforce := zforce - sum(ee.facet ff where not fixed,
		  ff.tension*ee.length*sqrt(ff.x^2+ff.y^2)/area);
             };
             printf "zforce: %18.15f  (geometric)\n", zforce;
           }