#include <stdio.h>
-#define NCL 11
+#define NCL 9
+#define Ndf 4
#define DYN_ON_SPEED 33 /*20 knots*/
extern COCKPIT cockpit_;
-SCALAR interp(SCALAR *y_table, SCALAR *x_table, int Ntable, SCALAR x)
+ SCALAR CLadot;
+ SCALAR CLq;
+ SCALAR CLde;
+ SCALAR CLob;
+
+
+ SCALAR Cdob;
+ SCALAR Cda; /*Not used*/
+ SCALAR Cdde;
+
+ SCALAR Cma;
+ SCALAR Cmadot;
+ SCALAR Cmq;
+ SCALAR Cmob;
+ SCALAR Cmde;
+
+ SCALAR Clbeta;
+ SCALAR Clp;
+ SCALAR Clr;
+ SCALAR Clda;
+ SCALAR Cldr;
+
+ SCALAR Cnbeta;
+ SCALAR Cnp;
+ SCALAR Cnr;
+ SCALAR Cnda;
+ SCALAR Cndr;
+
+ SCALAR Cybeta;
+ SCALAR Cyp;
+ SCALAR Cyr;
+ SCALAR Cyda;
+ SCALAR Cydr;
+
+ /*nondimensionalization quantities*/
+ /*units here are ft and lbs */
+ SCALAR cbar; /*mean aero chord ft*/
+ SCALAR b; /*wing span ft */
+ SCALAR Sw; /*wing planform surface area ft^2*/
+ SCALAR rPiARe; /*reciprocal of Pi*AR*e*/
+ SCALAR lbare; /*elevator moment arm MAC*/
+
+ SCALAR Weight; /*lbs*/
+ SCALAR MaxTakeoffWeight,EmptyWeight;
+ SCALAR Cg; /*%MAC*/
+ SCALAR Zcg; /*%MAC*/
+
+
+ SCALAR CLwbh,CL,cm,cd,cn,cy,croll,cbar_2V,b_2V,qS,qScbar,qSb;
+ SCALAR CLo,Cdo,Cmo;
+
+ SCALAR F_X_wind,F_Y_wind,F_Z_wind;
+
+ SCALAR long_trim;
+
+
+ SCALAR elevator, aileron, rudder;
+
+
+ SCALAR Flap_Position;
+
+ int Flaps_In_Transit;
+
+static SCALAR interp(SCALAR *y_table, SCALAR *x_table, int Ntable, SCALAR x)
{
SCALAR slope;
int i=1;
}
-
-void aero( SCALAR dt, int Initialize ) {
+void c172_aero( SCALAR dt, int Initialize ) {
- static int init = 0;
-
+ // static int init = 0;
+ static int fi=0;
+ static SCALAR lastFlapHandle=0;
+ static SCALAR Ai;
static SCALAR trim_inc = 0.0002;
- static SCALAR alpha_ind[NCL]={-0.087,0,0.175,0.209,0.24,0.262,0.278,0.303,0.314,0.332,0.367};
- static SCALAR CLtable[NCL]={-0.14,0.31,1.21,1.376,1.51249,1.591,1.63,1.60878,1.53712,1.376,1.142};
+ static SCALAR alpha_ind[NCL]={-0.087,0,0.14,0.21,0.24,0.26,0.28,0.31,0.35};
+ static SCALAR CLtable[NCL]={-0.22,0.25,1.02,1.252,1.354,1.44,1.466,1.298,0.97};
+
+ static SCALAR flap_ind[Ndf]={0,10,20,30};
+ static SCALAR flap_times[Ndf]={0,4,2,2};
+ static SCALAR dCLf[Ndf]={0,0.20,0.30,0.35};
+ static SCALAR dCdf[Ndf]={0,0.0021,0.0085,0.0191};
+ static SCALAR dCmf[Ndf]={0,-0.0654,-0.0981,-0.114};
+
+ static SCALAR flap_transit_rate=2.5;
+
+
+
/* printf("Initialize= %d\n",Initialize); */
/* printf("Initializing aero model...Initialize= %d\n", Initialize);
- */ CLadot=1.7;
+ */
+ /*nondimensionalization quantities*/
+ /*units here are ft and lbs */
+ cbar=4.9; /*mean aero chord ft*/
+ b=35.8; /*wing span ft */
+ Sw=174; /*wing planform surface area ft^2*/
+ rPiARe=0.054; /*reciprocal of Pi*AR*e*/
+ lbare=3.682; /*elevator moment arm / MAC*/
+
+ CLadot=1.7;
CLq=3.9;
- CLde=0.43;
- CLo=0;
-
+
+ CLob=0;
- Cdo=0.031;
+ Ai=1.24;
+ Cdob=0.036;
Cda=0.13; /*Not used*/
Cdde=0.06;
- Cma=-0.89;
+ Cma=-1.8;
Cmadot=-5.2;
Cmq=-12.4;
- Cmo=-0.062;
+ Cmob=-0.02;
Cmde=-1.28;
+
+ CLde=-Cmde / lbare; /* kinda backwards, huh? */
Clbeta=-0.089;
Clp=-0.47;
Clr=0.096;
- Clda=-0.178;
+ Clda=-0.09;
Cldr=0.0147;
Cnbeta=0.065;
Cnp=-0.03;
Cnr=-0.099;
- Cnda=-0.053;
+ Cnda=-0.0053;
Cndr=-0.0657;
Cybeta=-0.31;
Cyda=0.0;
Cydr=0.187;
- /*nondimensionalization quantities*/
- /*units here are ft and lbs */
- cbar=4.9; /*mean aero chord ft*/
- b=35.8; /*wing span ft */
- Sw=174; /*wing planform surface area ft^2*/
- rPiARe=0.054; /*reciprocal of Pi*AR*e*/
+
MaxTakeoffWeight=2550;
EmptyWeight=1500;
{ Weight=1500; }
- if(Dx_cg > 0.5586)
- { Dx_cg = 0.5586; }
- else if(Dx_cg < -0.4655)
- { Dx_cg = -0.4655; }
+ if(Dx_cg > 0.43)
+ { Dx_cg = 0.43; }
+ else if(Dx_cg < -0.6)
+ { Dx_cg = -0.6; }
- Cg=Dx_cg/cbar +0.25;
+ Cg=0.25 - Dx_cg/cbar;
Dz_cg=Zcg*cbar;
+ Dy_cg=0;
-
-
-
+
+ if(Flap_handle < flap_ind[0])
+ {
+ fi=0;
+ Flap_handle=flap_ind[0];
+ lastFlapHandle=Flap_handle;
+ Flap_Position=flap_ind[0];
+ }
+ else if(Flap_handle > flap_ind[Ndf-1])
+ {
+ fi=Ndf-1;
+ Flap_handle=flap_ind[fi];
+ lastFlapHandle=Flap_handle;
+ Flap_Position=flap_ind[fi];
+ }
+ else
+ {
+
+ if(dt <= 0)
+ Flap_Position=Flap_handle;
+ else
+ {
+ if(Flap_handle != lastFlapHandle)
+ {
+ Flaps_In_Transit=1;
+ }
+ if(Flaps_In_Transit)
+ {
+ fi=0;
+ while(flap_ind[fi] < Flap_handle) { fi++; }
+ if(Flap_Position < Flap_handle)
+ {
+ if(flap_times[fi] > 0)
+ flap_transit_rate=(flap_ind[fi] - flap_ind[fi-1])/flap_times[fi];
+ else
+ flap_transit_rate=(flap_ind[fi] - flap_ind[fi-1])/5;
+ }
+ else
+ {
+ if(flap_times[fi+1] > 0)
+ flap_transit_rate=(flap_ind[fi] - flap_ind[fi+1])/flap_times[fi+1];
+ else
+ flap_transit_rate=(flap_ind[fi] - flap_ind[fi+1])/5;
+ }
+ if(fabs(Flap_Position - Flap_handle) > dt*flap_transit_rate)
+ Flap_Position+=flap_transit_rate*dt;
+ else
+ {
+ Flaps_In_Transit=0;
+ Flap_Position=Flap_handle;
+ }
+ }
+ }
+ lastFlapHandle=Flap_handle;
+ }
- long_trim=0;
if(Aft_trim) long_trim = long_trim - trim_inc;
if(Fwd_trim) long_trim = long_trim + trim_inc;
/* printf("Long_control: %7.4f, long_trim: %7.4f,DEG_TO_RAD: %7.4f, RAD_TO_DEG: %7.4f\n",Long_control,long_trim,DEG_TO_RAD,RAD_TO_DEG);
*/ /*scale pct control to degrees deflection*/
- if ((Long_control+long_trim) <= 0)
- elevator=(Long_control+long_trim)*28*DEG_TO_RAD;
+ if ((Long_control+Long_trim) <= 0)
+ elevator=(Long_control+Long_trim)*28*DEG_TO_RAD;
else
- elevator=(Long_control+long_trim)*23*DEG_TO_RAD;
+ elevator=(Long_control+Long_trim)*23*DEG_TO_RAD;
aileron = -1*Lat_control*17.5*DEG_TO_RAD;
rudder = -1*Rudder_pedal*16*DEG_TO_RAD;
/*calculate rate derivative nondimensionalization (is that a word?) factors */
/*hack to avoid divide by zero*/
- /*the dynamic terms might be negligible at low ground speeds anyway*/
+ /*the dynamic terms are negligible at low ground speeds anyway*/
+/* printf("Vinf: %g, Span: %g\n",V_rel_wind,b);
+ */
if(V_rel_wind > DYN_ON_SPEED)
{
cbar_2V=cbar/(2*V_rel_wind);
/* printf("aero: Wb: %7.4f, Ub: %7.4f, Alpha: %7.4f, elev: %7.4f, ail: %7.4f, rud: %7.4f, long_trim: %7.4f\n",W_body,U_body,Alpha*RAD_TO_DEG,elevator*RAD_TO_DEG,aileron*RAD_TO_DEG,rudder*RAD_TO_DEG,long_trim*RAD_TO_DEG);
- */ //printf("Theta: %7.4f, Gamma: %7.4f, Beta: %7.4f, Phi: %7.4f, Psi: %7.4f\n",Theta*RAD_TO_DEG,Gamma_vert_rad*RAD_TO_DEG,Beta*RAD_TO_DEG,Phi*RAD_TO_DEG,Psi*RAD_TO_DEG);
+ printf("aero: Theta: %7.4f, Gamma: %7.4f, Beta: %7.4f, Phi: %7.4f, Psi: %7.4f\n",Theta*RAD_TO_DEG,Gamma_vert_rad*RAD_TO_DEG,Beta*RAD_TO_DEG,Phi*RAD_TO_DEG,Psi*RAD_TO_DEG);
+ */
+
/* sum coefficients */
- CLwbh = interp(CLtable,alpha_ind,NCL,Alpha);
- CL = CLo + CLwbh + (CLadot*Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
- cd = Cdo + rPiARe*CL*CL + Cdde*elevator;
- cy = Cybeta*Beta + (Cyp*P_body + Cyr*R_body)*b_2V + Cyda*aileron + Cydr*rudder;
+ CLwbh = interp(CLtable,alpha_ind,NCL,Std_Alpha);
+/* printf("CLwbh: %g\n",CLwbh);
+ */
+ CLo = CLob + interp(dCLf,flap_ind,Ndf,Flap_Position);
+ Cdo = Cdob + interp(dCdf,flap_ind,Ndf,Flap_Position);
+ Cmo = Cmob + interp(dCmf,flap_ind,Ndf,Flap_Position);
- cm = Cmo + Cma*Alpha + (Cmq*Q_body + Cmadot*Alpha_dot)*cbar_2V + Cmde*(elevator+long_trim);
- cn = Cnbeta*Beta + (Cnp*P_body + Cnr*R_body)*b_2V + Cnda*aileron + Cndr*rudder;
- croll=Clbeta*Beta + (Clp*P_body + Clr*R_body)*b_2V + Clda*aileron + Cldr*rudder;
+ /* printf("FP: %g\n",Flap_Position);
+ printf("CLo: %g\n",CLo);
+ printf("Cdo: %g\n",Cdo);
+ printf("Cmo: %g\n",Cmo); */
+
+
+
+
+
+ CL = CLo + CLwbh + (CLadot*Std_Alpha_dot + CLq*Theta_dot)*cbar_2V + CLde*elevator;
+ cd = Cdo + rPiARe*Ai*Ai*CL*CL + Cdde*elevator;
+ cy = Cybeta*Std_Beta + (Cyp*P_body + Cyr*R_body)*b_2V + Cyda*aileron + Cydr*rudder;
+
+ cm = Cmo + Cma*Std_Alpha + (Cmq*Q_body + Cmadot*Std_Alpha_dot)*cbar_2V + Cmde*(elevator);
+ cn = Cnbeta*Std_Beta + (Cnp*P_body + Cnr*R_body)*b_2V + Cnda*aileron + Cndr*rudder;
+ croll=Clbeta*Std_Beta + (Clp*P_body + Clr*R_body)*b_2V + Clda*aileron + Cldr*rudder;
/* printf("aero: CL: %7.4f, Cd: %7.4f, Cm: %7.4f, Cy: %7.4f, Cn: %7.4f, Cl: %7.4f\n",CL,cd,cm,cy,cn,croll);
- */ /*calculate wind axes forces*/
+ */
+
+ /*calculate wind axes forces*/
F_X_wind=-1*cd*qS;
F_Y_wind=cy*qS;
F_Z_wind=-1*CL*qS;
/* printf("V_rel_wind: %7.4f, Fxwind: %7.4f Fywind: %7.4f Fzwind: %7.4f\n",V_rel_wind,F_X_wind,F_Y_wind,F_Z_wind);
- */
+ */
+
/*calculate moments and body axis forces */
M_n_aero = cn*qSb;
/* printf("I_yy: %7.4f, qScbar: %7.4f, qbar: %7.4f, Sw: %7.4f, cbar: %7.4f, 0.5*rho*V^2: %7.4f\n",I_yy,qScbar,Dynamic_pressure,Sw,cbar,0.5*0.0023081*V_rel_wind*V_rel_wind);
- */
-/* printf("Fxaero: %7.4f Fyaero: %7.4f Fzaero: %7.4f Weight: %7.4f\n",F_X_aero,F_Y_aero,F_Z_aero,W);
- *//* printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);
- */
+
+ printf("Fxaero: %7.4f Fyaero: %7.4f Fzaero: %7.4f Weight: %7.4f\n",F_X_aero,F_Y_aero,F_Z_aero,Weight);
+
+ printf("Maero: %7.4f Naero: %7.4f Raero: %7.4f\n",M_m_aero,M_n_aero,M_l_aero);
+ */
+
}
projects / flightgear.git / blobdiff