// Commands to calculate x, y, and z forces, and tilt torque. // command to smoothly change x offset. // Use: set new_x_offset to desired value, then do do_x_offset. new_x_offset := x_offset; do_x_offset := { doff := new_x_offset - x_offset; x_offset := new_x_offset; slope := tan(tilt*pi/180); set vertex x x + doff*z/(height + slope*(y - y_offset)); recalc; } // Central difference, moving linearly dx := 1e-5; do_xforce := { new_x_offset := x_offset + dx; // do this first so constraint in right place do_x_offset; new_energy_1 := total_energy - body[1].pressure * (body[1].volume - body[1].target); new_x_offset := x_offset - 2*dx; do_x_offset; new_energy_2 := total_energy - body[1].pressure * (body[1].volume - body[1].target); xforce := -(new_energy_1 - new_energy_2)/2/dx; printf "xforce: %18.15f (central difference, linear move)\n", xforce; new_x_offset := x_offset + dx; // restore original do_x_offset; } // command to smoothly change y offset. // Use: set new_y_offset to desired value, then do do_y_offset. new_y_offset := y_offset; do_y_offset := { doff := new_y_offset - y_offset; y_offset := new_y_offset; set vertex y y + doff*z/(height + slope*(y - y_offset)); recalc; } // Central difference, moving linearly dy := 1e-5; do_yforce := { new_y_offset := y_offset + dy; // do this first so constraint in right place do_y_offset; new_energy_1 := total_energy - body[1].pressure * (body[1].volume - body[1].target); new_y_offset := y_offset - 2*dy; do_y_offset; new_energy_2 := total_energy - body[1].pressure * (body[1].volume - body[1].target); yforce := -(new_energy_1 - new_energy_2)/2/dy; printf "yforce: %18.15f (central difference, linear move)\n", yforce; new_y_offset := y_offset + dy; // restore original do_y_offset; } // Central difference, moving linearly dz := 1e-5; do_zforce := { 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 + slope*(y - y_offset)); 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 + slope*(y - y_offset)); height := height + dz; // restore original recalc; } // Torque // First, a smooth tilting new_tilt := tilt; do_tilt := { dtilt := (new_tilt - tilt)*pi/180; slope := tan(tilt*pi/180); tilt := new_tilt; foreach vertex vv do { beta := z/(height + slope*(y - y_offset)); vv.y := y_offset + cos(beta*dtilt)*(y - y_offset) - sin(beta*dtilt)*(z - beta*height); vv.z := beta*height + sin(beta*dtilt)*(y - y_offset) + cos(beta*dtilt)*(z - beta*height); } } // Torque by central difference, moving linearly dangle := 1e-4; // degrees do_torque := { new_tilt := tilt + dangle; // do this first so constraint in right place do_tilt; new_energy_1 := total_energy - body[1].pressure * (body[1].volume - body[1].target); new_tilt := tilt - 2*dangle; do_tilt; new_energy_2 := total_energy - body[1].pressure * (body[1].volume - body[1].target); torque := -(new_energy_1 - new_energy_2)/2/(dangle*pi/180); printf "torque: %18.15f (central difference, linear move)\n", torque; new_tilt := tilt + dangle; // restore original do_tilt; } forces := { do_xforce; do_yforce; do_zforce; do_torque; }