Index: mbdyn/aero/genfm.cc
===================================================================
RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/genfm.cc,v
retrieving revision 1.11
retrieving revision 1.15
diff -u -r1.11 -r1.15
--- mbdyn/aero/genfm.cc	13 Feb 2010 00:00:26 -0000	1.11
+++ mbdyn/aero/genfm.cc	21 Jul 2010 13:55:00 -0000	1.15
@@ -33,6 +33,8 @@
 #include "mbconfig.h"           /* This goes first in every *.c,*.cc file */
 #endif /* HAVE_CONFIG_H */
 
+#include <cmath>
+
 #include "dataman.h"
 #include "aeroelem.h"
 #include "genfm.h"
@@ -68,19 +70,17 @@
 	}
 #endif
 
-	Vec3 V(R.Transpose()*Vca);
+	Vec3 V(R.MulTV(Vca));
 
-	/* FIXME: according to source, the rotation sequence
-	 * is alpha, then beta */
-	dAlpha = atan2(V(3), V(1));
-	dBeta = -atan2(V(2), V(1));
+	dAlpha = atan2(V(3), std::abs(V(1)));
+	dBeta = atan2(-V(2), copysign(std::sqrt(V(1)*V(1) + V(2)*V(2)), V(1)));
 
 	doublereal rho, c, p, T;
 	GetAirProps(Xca, rho, c, p, T);	/* p, T no used yet */
 	doublereal q = Vca.Dot()*rho/2.;
 
 	doublereal dScaleForce = q*dRefSurface;
-	doublereal dScaleMoment = dScaleForce*dRefSurface;
+	doublereal dScaleMoment = dScaleForce*dRefLength;
 
 	int nAlpha = pData->Alpha.size() - 1;
 	int nBeta = pData->Beta.size() - 1;
@@ -88,27 +88,26 @@
 	ASSERT(nAlpha >= 0);
 	ASSERT(nBeta >= 0);
 
-	/* TODO: search for Alpha, Beta */
-	if (dAlpha <= pData->Alpha[0]) {
+	if (dAlpha < pData->Alpha[0]) {
 		/* smooth out coefficients if Alpha does not span -180 => 180 */
 		doublereal dAlphaX = (dAlpha - pData->Alpha[0])/(-M_PI - pData->Alpha[0]);
 		doublereal dSmoothAlpha = (std::cos(M_PI*dAlphaX) + 1)/2.;
 
-		if (dBeta <= pData->Beta[0]) {
+		if (dBeta < pData->Beta[0]) {
 			/* smooth out coefficients if Beta does not span -180 => 180 */
 			doublereal dBetaX = (dBeta - pData->Beta[0])/(-M_PI - pData->Beta[0]);
 			doublereal dSmoothBeta = (std::cos(M_PI*dBetaX) + 1)/2.;
 
-			F = Vec3(&pData->Data[0][0].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
-			M = Vec3(&pData->Data[0][0].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
+			tilde_F = Vec3(&pData->Data[0][0].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
+			tilde_M = Vec3(&pData->Data[0][0].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
 
-		} else if (dBeta >= pData->Beta[nBeta]) {
+		} else if (dBeta > pData->Beta[nBeta]) {
 			/* smooth out coefficients if Beta does not span -180 => 180 */
 			doublereal dBetaX = (dBeta - pData->Beta[nBeta])/(M_PI - pData->Beta[nBeta]);
 			doublereal dSmoothBeta = (std::cos(M_PI*dBetaX) + 1)/2.;
 
-			F = Vec3(&pData->Data[nBeta][0].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
-			M = Vec3(&pData->Data[nBeta][0].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
+			tilde_F = Vec3(&pData->Data[nBeta][0].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
+			tilde_M = Vec3(&pData->Data[nBeta][0].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
 
 		} else {
 			int iBeta = bisec<doublereal>(&pData->Beta[0], dBeta, 0, nBeta);
@@ -123,30 +122,30 @@
 			GenericAerodynamicData::GenericAerodynamicCoef c
 				= pData->Data[iBeta][0]*d1Beta + pData->Data[iBeta + 1][0]*d2Beta;
 
-			F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothAlpha);
-			M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothAlpha);
+			tilde_F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothAlpha);
+			tilde_M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothAlpha);
 		}
 
-	} else if (dAlpha >= pData->Alpha[nAlpha]) {
+	} else if (dAlpha > pData->Alpha[nAlpha]) {
 		/* smooth out coefficients if Alpha does not span -180 => 180 */
 		doublereal dAlphaX = (dAlpha - pData->Alpha[nAlpha])/(-M_PI - pData->Alpha[nAlpha]);
 		doublereal dSmoothAlpha = (std::cos(M_PI*dAlphaX) + 1)/2.;
 
-		if (dBeta <= pData->Beta[0]) {
+		if (dBeta < pData->Beta[0]) {
 			/* smooth out coefficients if Beta does not span -180 => 180 */
 			doublereal dBetaX = (dBeta - pData->Beta[0])/(-M_PI - pData->Beta[0]);
 			doublereal dSmoothBeta = (std::cos(M_PI*dBetaX) + 1)/2.;
 
-			F = Vec3(&pData->Data[0][nAlpha].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
-			M = Vec3(&pData->Data[0][nAlpha].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
+			tilde_F = Vec3(&pData->Data[0][nAlpha].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
+			tilde_M = Vec3(&pData->Data[0][nAlpha].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
 
-		} else if (dBeta >= pData->Beta[nBeta]) {
+		} else if (dBeta > pData->Beta[nBeta]) {
 			/* smooth out coefficients if Beta does not span -180 => 180 */
 			doublereal dBetaX = (dBeta - pData->Beta[nBeta])/(M_PI - pData->Beta[nBeta]);
 			doublereal dSmoothBeta = (std::cos(M_PI*dBetaX) + 1)/2.;
 
-			F = Vec3(&pData->Data[nBeta][nAlpha].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
-			M = Vec3(&pData->Data[nBeta][nAlpha].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
+			tilde_F = Vec3(&pData->Data[nBeta][nAlpha].dCoef[0])*(dScaleForce*dSmoothAlpha*dSmoothBeta);
+			tilde_M = Vec3(&pData->Data[nBeta][nAlpha].dCoef[3])*(dScaleMoment*dSmoothAlpha*dSmoothBeta);
 
 		} else {
 			int iBeta = bisec<doublereal>(&pData->Beta[0], dBeta, 0, nBeta);
@@ -161,8 +160,8 @@
 			GenericAerodynamicData::GenericAerodynamicCoef c
 				= pData->Data[iBeta][nAlpha]*d1Beta + pData->Data[iBeta + 1][nAlpha]*d2Beta;
 
-			F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothAlpha);
-			M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothAlpha);
+			tilde_F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothAlpha);
+			tilde_M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothAlpha);
 		}
 
 	} else {
@@ -175,7 +174,7 @@
 		doublereal d1Alpha = (pData->Alpha[iAlpha + 1] - dAlpha)/ddAlpha;
 		doublereal d2Alpha = (dAlpha - pData->Alpha[iAlpha])/ddAlpha;
 
-		if (dBeta <= pData->Beta[0]) {
+		if (dBeta < pData->Beta[0]) {
 			/* smooth out coefficients if Beta does not span -180 => 180 */
 			doublereal dBetaX = (dBeta - pData->Beta[0])/(-M_PI - pData->Beta[0]);
 			doublereal dSmoothBeta = (std::cos(M_PI*dBetaX) + 1)/2.;
@@ -183,10 +182,10 @@
 			GenericAerodynamicData::GenericAerodynamicCoef c
 				= pData->Data[0][iAlpha]*d1Alpha + pData->Data[0][iAlpha + 1]*d2Alpha;
 
-			F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothBeta);
-			M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothBeta);
+			tilde_F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothBeta);
+			tilde_M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothBeta);
 
-		} else if (dBeta >= pData->Beta[nBeta]) {
+		} else if (dBeta > pData->Beta[nBeta]) {
 			/* smooth out coefficients if Beta does not span -180 => 180 */
 			doublereal dBetaX = (dBeta - pData->Beta[nBeta])/(M_PI - pData->Beta[nBeta]);
 			doublereal dSmoothBeta = (std::cos(M_PI*dBetaX) + 1)/2.;
@@ -194,8 +193,8 @@
 			GenericAerodynamicData::GenericAerodynamicCoef c
 				= pData->Data[nBeta][iAlpha]*d1Alpha + pData->Data[nBeta][iAlpha + 1]*d2Alpha;
 
-			F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothBeta);
-			M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothBeta);
+			tilde_F = Vec3(&c.dCoef[0])*(dScaleForce*dSmoothBeta);
+			tilde_M = Vec3(&c.dCoef[3])*(dScaleMoment*dSmoothBeta);
 
 		} else {
 			int iBeta = bisec<doublereal>(&pData->Beta[0], dBeta, 0, nBeta);
@@ -214,11 +213,14 @@
 			
 			GenericAerodynamicData::GenericAerodynamicCoef c = c1*d1Beta + c2*d2Beta;
 
-			F = Vec3(&c.dCoef[0])*dScaleForce;
-			M = Vec3(&c.dCoef[3])*dScaleMoment;
+			tilde_F = Vec3(&c.dCoef[0])*dScaleForce;
+			tilde_M = Vec3(&c.dCoef[3])*dScaleMoment;
 		}
 	}
 
+	F = R*tilde_F;
+	M = R*tilde_M + f.Cross(F);
+
 	WorkVec.Add(1, F);
 	WorkVec.Add(4, M);
 }
@@ -238,6 +240,8 @@
 dRefLength(dL),
 tilde_f(fTmp),
 tilde_Ra(RaTmp),
+tilde_F(0.),
+tilde_M(0.),
 F(0.),
 M(0.),
 pData(pD)
@@ -252,6 +256,72 @@
 	}
 }
 
+/* Tipo dell'elemento (usato per debug ecc.) */
+Elem::Type
+GenericAerodynamicForce::GetElemType(void) const
+{
+	return Elem::AERODYNAMIC;
+}
+
+/* Dimensioni del workspace */
+void
+GenericAerodynamicForce::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
+{
+	*piNumRows = 6;
+	*piNumCols = 1;
+}
+
+/* assemblaggio jacobiano */
+VariableSubMatrixHandler&
+GenericAerodynamicForce::AssJac(VariableSubMatrixHandler& WorkMat,
+	doublereal /* dCoef */ ,
+	const VectorHandler& /* XCurr */ ,
+	const VectorHandler& /* XPrimeCurr */ )
+{
+	DEBUGCOUTFNAME("GenericAerodynamicForce::AssJac");
+	WorkMat.SetNullMatrix();
+	return WorkMat;
+}
+
+/* Numero di GDL iniziali */
+unsigned int
+GenericAerodynamicForce::iGetInitialNumDof(void) const
+{
+	return 0;
+}
+
+/* Dimensioni del workspace */
+void
+GenericAerodynamicForce::InitialWorkSpaceDim(integer* piNumRows, integer* piNumCols) const
+{
+	*piNumRows = 6;
+	*piNumCols = 1;
+}
+
+/* assemblaggio jacobiano */
+VariableSubMatrixHandler&
+GenericAerodynamicForce::InitialAssJac(VariableSubMatrixHandler& WorkMat,
+	const VectorHandler& /* XCurr */)
+{
+	DEBUGCOUTFNAME("GenericAerodynamicForce::InitialAssJac");
+	WorkMat.SetNullMatrix();
+	return WorkMat;
+}
+
+/* Tipo di elemento aerodinamico */
+AerodynamicElem::Type
+GenericAerodynamicForce::GetAerodynamicElemType(void) const
+{
+	return AerodynamicElem::GENERICFORCE;
+}
+
+void
+GenericAerodynamicForce::GetConnectedNodes(std::vector<const Node *>& connectedNodes) const
+{
+	connectedNodes.resize(1);
+	connectedNodes[0] = pNode;
+}
+
 /* Scrive il contributo dell'elemento al file di restart */
 std::ostream&
 GenericAerodynamicForce::Restart(std::ostream& out) const
@@ -290,6 +360,7 @@
 			<< std::setw(8) << GetLabel()
 			<< " " << dAlpha*180./M_PI
 			<< " " << dBeta*180./M_PI
+			<< " " << tilde_F << " " << tilde_M
 			<< " " << F << " " << M << std::endl;
 	}
 }
@@ -581,11 +652,17 @@
 
 	/* The offset is in the reference frame of the node */
 	ReferenceFrame RF(pNode);
-	Vec3 f(HP.GetPosRel(RF));
+	Vec3 f(0.);
+	if (HP.IsKeyWord("position")) {
+		f = HP.GetPosRel(RF);
+	}
 
 	/* the orientation is in flight mechanics (FIXME?) reference frame:
 	 * X forward, Y to the right, Z down */
-	Mat3x3 Ra(HP.GetRotRel(RF));
+	Mat3x3 Ra(Eye3);
+	if (HP.IsKeyWord("orientation")) {
+		Ra = HP.GetRotRel(RF);
+	}
 
 	/* 1. by default, which means that coefficients are only normalized
 	 * by the dynamic pressure */
@@ -601,8 +678,23 @@
 	}
 
 	/* TODO: allow to reference previously loaded data */
-	std::string fname(HP.GetFileName());
-	GenericAerodynamicData *pData = ReadGenericAerodynamicData(fname);
+	GenericAerodynamicData *pData = 0;
+	if (HP.IsKeyWord("file")) {
+		std::string fname(HP.GetFileName());
+		pData = ReadGenericAerodynamicData(fname);
+
+	} else if (HP.IsKeyWord("reference")) {
+		silent_cerr("GenericAerodynamicForce(" << uLabel << "): "
+			"references to generic aerodynamic data not implemented yet "
+			"at line " << HP.GetLineData() << std::endl);
+		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
+
+	} else {
+		silent_cerr("GenericAerodynamicForce(" << uLabel << "): "
+			"keyword \"file\" expected "
+			"at line " << HP.GetLineData() << std::endl);
+		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
+	}
 
 	flag fOut = pDM->fReadOutput(HP, Elem::AERODYNAMIC);
 
Index: mbdyn/aero/genfm.h
===================================================================
RCS file: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/genfm.h,v
retrieving revision 1.8
retrieving revision 1.10
diff -u -r1.8 -r1.10
--- mbdyn/aero/genfm.h	13 Feb 2010 00:00:26 -0000	1.8
+++ mbdyn/aero/genfm.h	21 Jul 2010 13:55:00 -0000	1.10
@@ -82,6 +82,8 @@
 	const Mat3x3 tilde_Ra;	/* Rotaz. del sistema aerodinamico al nodo */
 
 	/* force and moment */
+	Vec3 tilde_F;
+	Vec3 tilde_M;
 	Vec3 F;
 	Vec3 M;
 
@@ -107,31 +109,20 @@
 	virtual std::ostream& Restart(std::ostream& out) const;
 
 	/* Tipo dell'elemento (usato per debug ecc.) */
-	virtual Elem::Type GetElemType(void) const {
-		return Elem::AERODYNAMIC;
-	};
+	virtual Elem::Type GetElemType(void) const;
 
 	/* funzioni proprie */
 
 	/* Dimensioni del workspace */
 	virtual void
-	WorkSpaceDim(integer* piNumRows, integer* piNumCols) const {
-		*piNumRows = 6;
-		*piNumCols = 1;
-		/* NOTE: Jacobian could be implemented... */
-	};
+	WorkSpaceDim(integer* piNumRows, integer* piNumCols) const;
 
 	/* assemblaggio jacobiano */
 	virtual VariableSubMatrixHandler&
 	AssJac(VariableSubMatrixHandler& WorkMat,
-		doublereal /* dCoef */ ,
-		const VectorHandler& /* XCurr */ ,
-		const VectorHandler& /* XPrimeCurr */ )
-	{
-		DEBUGCOUTFNAME("GenericAerodynamicForce::AssJac");
-		WorkMat.SetNullMatrix();
-		return WorkMat;
-	};
+		doublereal dCoef,
+		const VectorHandler& XCurr,
+		const VectorHandler& XPrimeCurr);
 
 	/* assemblaggio residuo */
 	virtual SubVectorHandler&
@@ -147,28 +138,16 @@
 	virtual void Output(OutputHandler& OH) const;
 
 	/* Numero di GDL iniziali */
-	virtual unsigned int iGetInitialNumDof(void) const
-	{
-		return 0;
-	};
+	virtual unsigned int iGetInitialNumDof(void) const;
 
 	/* Dimensioni del workspace */
 	virtual void
-	InitialWorkSpaceDim(integer* piNumRows, integer* piNumCols) const
-	{
-		*piNumRows = 6;
-		*piNumCols = 1;
-	};
+	InitialWorkSpaceDim(integer* piNumRows, integer* piNumCols) const;
 
 	/* assemblaggio jacobiano */
 	virtual VariableSubMatrixHandler&
 	InitialAssJac(VariableSubMatrixHandler& WorkMat,
-		const VectorHandler& /* XCurr */)
-	{
-		DEBUGCOUTFNAME("GenericAerodynamicForce::InitialAssJac");
-		WorkMat.SetNullMatrix();
-		return WorkMat;
-	};
+		const VectorHandler& XCurr);
 
 	/* assemblaggio residuo */
 	virtual SubVectorHandler&
@@ -176,22 +155,14 @@
 		const VectorHandler& XCurr);
 
 	/* Tipo di elemento aerodinamico */
-	virtual AerodynamicElem::Type GetAerodynamicElemType(void) const
-	{
-		return AerodynamicElem::GENERICFORCE;
-	};
+	virtual AerodynamicElem::Type GetAerodynamicElemType(void) const;
 
-	/* *******PER IL SOLUTORE PARALLELO******** */
 	/*
 	 * Fornisce il tipo e la label dei nodi che sono connessi all'elemento
 	 * utile per l'assemblaggio della matrice di connessione fra i dofs
 	 */
 	virtual void
-	GetConnectedNodes(std::vector<const Node *>& connectedNodes) const {
-		connectedNodes.resize(1);
-		connectedNodes[0] = pNode;
-	};
-	/* ************************************************ */
+	GetConnectedNodes(std::vector<const Node *>& connectedNodes) const;
 };
 
 extern Elem *