// gullforces.cmd

// Commands to calculate forces on the gull wing lead.

//***************************************************************************

// Auxiliary command to calculate relative distance from pad and lead
// for purposes of interpolating movement.  The vertex extra parameter alpha
// is set to 1 on the lead, 0 on the pad, and linearly interpolate in
// z up to the lowest contact line vertex on the lead.
// Things get so tricky in order to handle moving a lead that has gap = 0.

define vertex attribute alpha real

calc_alpha := { 
   minz := min(vertex where on_constraint 13 or on_constraint 23, z);
   if ( minz <= 0.0 ) then
      printf "WARNING: contact line hits pad.\n";
   foreach vertex vv do
   { if ( vv.on_constraint 1 or vv.on_constraint 5 ) then
       { vv.alpha := 0; continue; };
     if ( vv.fixed ) then { vv.alpha := 1; continue; };
     vv.alpha := vv.z >= minz ?  1 : vv.z/minz ;
   }
}

//************************************************************************
// Some lead-moving commands.  These use the alpha attribute to 
// interpolate motion.

// Move in x direction by delta_x, i.e. off-center 
delta_x := 0
move_x := { calc_alpha;
            leadshift := leadshift + delta_x;  // do parameters first
            set vertex x x+alpha*delta_x;
          }

// Move in y direction by delta_y
delta_y := 0
move_y := { calc_alpha;
            leadheel := leadheel + delta_y;
            leadtoe := leadtoe + delta_y; 
            set vertex y y+alpha*delta_y;
          }

// Move in z direction by delta_z
delta_z := 0
move_z := { calc_alpha;
            gap := gap + delta_z; 
            set vertex z z+alpha*delta_z;
          }

//**************************************************************************

// Force calculations by central differences.

dx := 0.0001
calc_xforce := { 
   delta_x := dx;
   move_x;
   hi := total_energy - body[1].pressure*(body[1].volume-body[1].target);
   delta_x := -2*dx;
   move_x;
   lo := total_energy - body[1].pressure*(body[1].volume-body[1].target);
   delta_x := dx;
   move_x;
   xforce := -(hi - lo)/2/dx;
}

dy := 0.0001
calc_yforce := { 
   delta_y := dy;
   move_y;
   hi := total_energy - body[1].pressure*(body[1].volume-body[1].target);
   delta_y := -2*dy;
   move_y;
   lo := total_energy - body[1].pressure*(body[1].volume-body[1].target);
   delta_y := dy;
   move_y;
   yforce := -(hi - lo)/2/dy;
}

dz := 0.0001
calc_zforce := { 
   delta_z := dz;
   move_z;
   hi := total_energy - body[1].pressure*(body[1].volume-body[1].target);
   delta_z := -2*dz;
   move_z;
   lo := total_energy - body[1].pressure*(body[1].volume-body[1].target);
   delta_z := dz;
   move_z;
   zforce := -(hi - lo)/2/dz;
}

forces := { calc_xforce; printf "xforce: %18.14f\n",xforce;
            calc_yforce; printf "yforce: %18.14f\n",yforce;
            calc_zforce; printf "zforce: %18.14f\n",zforce;
          }