RESSORTS

//-------------------------------------------------
typedef struct tagRESSORT
  {
  double    k;   // constante de raideur
  double    l0;  // longeur a vide
  }RESSORT,*P_RESSORT,**PP_RESSORT;

//  il y a 3 masses A, B, et C
// liees par des ressorts
typedef struct tagEQUATION
  {
  // les donnees
  double    mA,mB,mC;       // masses
  RESSORT   rAB,rAC,rBC;    // les trois ressorts  
  
  // informations pour le dessin  
  COLORREF  cA,cB,cC;       // couleurs des masses
  COLORREF  cRessorts;      // couleurs des ressorts
  DWORD     t0;

  double    alpha;          // facteur d'amortissement
  
  // l'integrateur d'equations differentielles
  ODE_RK  rk;
  }EQUATION,*P_EQUATION,**PP_EQUATION;


//-------------------------------------------------
// 12 valeurs : (3 masses) * (position (x,y) + vitesse (vx,vy)) = (3)*(2+2) = 12
//-------------------------------------------------
#define INDEX_XA        0
#define INDEX_YA        1
#define INDEX_VXA       2
#define INDEX_VYA       3

#define INDEX_XB        4
#define INDEX_YB        5
#define INDEX_VXB       6
#define INDEX_VYB       7

#define INDEX_XC        8
#define INDEX_YC        9
#define INDEX_VXC       10
#define INDEX_VYC       11
//-------------------------------------------------
// 12 valeurs : (3 masses) * (vitesse (vx,vy) + acceleration (ax,ay)) = (3)*(2+2) = 12
//-------------------------------------------------
#define DINDEX_VXA       0
#define DINDEX_VYA       1
#define DINDEX_AXA       2
#define DINDEX_AYA       3

#define DINDEX_VXB       4
#define DINDEX_VYB       5
#define DINDEX_AXB       6
#define DINDEX_AYB       7

#define DINDEX_VXC       8
#define DINDEX_VYC       9
#define DINDEX_AXC       10
#define DINDEX_AYC       11

//-------------------------------------------------
// on cherche a resoudre Y'=f(X,Y)
// avec X un scalaire, et Y et Y' sont des vecteurs (12 dimensions) :
// Y = ( xA  )  Y' = ( vxA )  <-- masse A
//     ( yA  )       ( vyA )  
//     ( vxA )       ( axA )  
//     ( vyA )       ( ayA )  
//     ( xB  )       ( vxB )  <-- masse B
//     ( yB  )       ( vyB )  
//     ( vxB )       ( axB )  
//     ( vyB )       ( ayB )  
//     ( xC  )       ( vxC )  <-- masse C
//     ( yC  )       ( vyC )  
//     ( vxC )       ( axC )  
//     ( vyC )       ( ayC )  
// 
// ici dans cette fonction <res> est le vecteur Y' que l'on doit calculer
BOOL funcEquaDiff(P_VECTND res,double x,P_VECTND y,void *param)
{
P_EQUATION  eq;
double      FAB,FAC,FBC;              // les forces (en norme relative +/-)
double      lAB,lAC,lBC;              // les longueurs
double      angleAB,angleAC,angleBC;  // les angles des droites (AB), (AC), et (BC)
double      dxAB,dyAB;                // variables utiles (calcul des longueurs et des angles)
double      dxAC,dyAC;                // variables utiles (calcul des longueurs et des angles)
double      dxBC,dyBC;                // variables utiles (calcul des longueurs et des angles)

eq = (P_EQUATION)param;
AssertPointer(eq);

// le probleme des vitesses est trivial,
// puisque l'on nous les donne
res->coord[DINDEX_VXA] = y->coord[INDEX_VXA]; // masse A
res->coord[DINDEX_VYA] = y->coord[INDEX_VYA];
res->coord[DINDEX_VXB] = y->coord[INDEX_VXB]; // masse B
res->coord[DINDEX_VYB] = y->coord[INDEX_VYB];
res->coord[DINDEX_VXC] = y->coord[INDEX_VXC]; // masse C
res->coord[DINDEX_VYC] = y->coord[INDEX_VYC];

// passons aux choses serieuses, l'accelerations des masses
// est donnees par le PFD (principe fondamental de la dynamique)
// a=F/m, de plus on sait calculer F pour des ressorts : F = - k.(l-l0)

// * on calcule les longueurs
dxAB = y->coord[INDEX_XB] - y->coord[INDEX_XA]; // longueur entre A et B
dyAB = y->coord[INDEX_YB] - y->coord[INDEX_YA];
lAB  = sqrt(dxAB*dxAB + dyAB*dyAB);
dxAC = y->coord[INDEX_XC] - y->coord[INDEX_XA]; // longueur entre A et C
dyAC = y->coord[INDEX_YC] - y->coord[INDEX_YA];
lAC  = sqrt(dxAC*dxAC + dyAC*dyAC);
dxBC = y->coord[INDEX_XC] - y->coord[INDEX_XB]; // longueur entre B et C
dyBC = y->coord[INDEX_YC] - y->coord[INDEX_YB];
lBC  = sqrt(dxBC*dxBC + dyBC*dyBC);
// * on calcule les forces
FAB = - eq->rAB.k * (lAB - eq->rAB.l0);
FAC = - eq->rAC.k * (lAC - eq->rAC.l0);
FBC = - eq->rBC.k * (lBC - eq->rBC.l0);
// * on calcule les angles (sens trigonometrique)
angleAB = atan2(dyAB,dxAB);
angleAC = atan2(dyAC,dxAC);
angleBC = atan2(dyBC,dxBC);
// * on calcule les accelerations (attention aux signes)
res->coord[DINDEX_AXA] = ((- FAB*cos(angleAB) - FAC*cos(angleAC))/eq->mA) - eq->alpha*y->coord[INDEX_VXA]; // masse A
res->coord[DINDEX_AYA] = ((- FAB*sin(angleAB) - FAC*sin(angleAC))/eq->mA) - eq->alpha*y->coord[INDEX_VYA];
res->coord[DINDEX_AXB] = ((+ FAB*cos(angleAB) - FBC*cos(angleBC))/eq->mB) - eq->alpha*y->coord[INDEX_VXB]; // masse B
res->coord[DINDEX_AYB] = ((+ FAB*sin(angleAB) - FBC*sin(angleBC))/eq->mB) - eq->alpha*y->coord[INDEX_VYB];
res->coord[DINDEX_AXC] = ((+ FAC*cos(angleAC) + FBC*cos(angleBC))/eq->mC) - eq->alpha*y->coord[INDEX_VXC]; // masse C
res->coord[DINDEX_AYC] = ((+ FAC*sin(angleAC) + FBC*sin(angleBC))/eq->mC) - eq->alpha*y->coord[INDEX_VYC];

return TRUE;
} // funcEquaDiff()