//
// Written by Curtis Olson, started October 1998.
//
-// Copyright (C) 1998 Curtis L. Olson - curt@me.umn.edu
+// Copyright (C) 1998 Curtis L. Olson - http://www.flightgear.org/~curt
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
//
// $Id$
#include <stdio.h>
#include <string.h>
-
-
-
void do_trims(int kmax,FILE *out,InitialConditions IC)
{
- int k=0,i,j;
+ int bad_trim=0,i,j;
double speed,elevator,cmcl,maxspeed;
out=fopen("trims.out","w");
speed=55;
- for(j=0;j<=30;j+=10)
+ for(j=0;j<=0;j+=10)
{
IC.flap_handle=j;
- for(i=1;i<=5;i++)
+ for(i=4;i<=4;i++)
{
switch(i)
{
case 1: IC.weight=1500;IC.cg=0.155;break;
case 2: IC.weight=1500;IC.cg=0.364;break;
case 3: IC.weight=1950;IC.cg=0.155;break;
- case 4: IC.weight=2550;IC.cg=0.257;break;
+ case 4: IC.weight=2400;IC.cg=0.257;break;
case 5: IC.weight=2550;IC.cg=0.364;break;
}
IC.vc=speed;
Long_control=0;Theta=0;Throttle_pct=0.0;
- k=trim_long(kmax,IC);
+ bad_trim=trim_long(kmax,IC);
if(Long_control <= 0)
elevator=Long_control*28;
else
{
cmcl=cm / CL;
}
- if(k < kmax)
+ if(!bad_trim)
{
- fprintf(out,"%g,%g,%g,%g,%g,%d",V_calibrated_kts,Alpha*RAD_TO_DEG,Long_control,Throttle_pct,Flap_Position,k);
+ fprintf(out,"%g,%g,%g,%g,%g",V_calibrated_kts,Alpha*RAD_TO_DEG,Long_control,Throttle_pct,Flap_Position);
fprintf(out,",%g,%g,%g,%g,%g\n",CL,cm,cmcl,Weight,Cg);
/* printf("%g,%g,%g,%g,%g,%g,%g,%g,%g,%g\n",V_calibrated_kts,Alpha*RAD_TO_DEG,elevator,CL,cm,Cmo,Cma,Cmde,Mass*32.174,Dx_cg);
*/ }
else
{
- printf("kmax exceeded at: %g knots, %g lbs, %g %%MAC, Flaps: %g\n",V_calibrated_kts,Weight,Cg,Flap_Position);
+ printf("kmax exceeded at: %g knots, %g lbs, %g %%MAC, Flaps: %g\n",V_true_kts,Weight,Cg,Flap_Position);
printf("wdot: %g, udot: %g, qdot: %g\n",W_dot_body,U_dot_body,Q_dot_body);
printf("Alpha: %g, Throttle_pct: %g, Long_control: %g\n\n",Alpha*RAD_TO_DEG,Throttle_pct,Long_control);
}
}
}
fclose(out);
-}
+}
+
+find_max_alt(int kmax,InitialConditions IC)
+{
+ int bad_trim=0,i=0;
+ float min=0,max=30000;
+ IC.use_gamma_tmg=1;
+ IC.gamma=0;
+ IC.vc=73;
+ IC.altitude==1000;
+ while(!bad_trim)
+ {
+ bad_trim=trim_long(200,IC);
+ IC.altitude+=1000;
+ }
+ while((fabs(max-min) > 100) && (i < 50))
+ {
+
+ IC.altitude=(max-min)/2 + min;
+ printf("\nIC.altitude: %g, max: %g, min: %g, bad_trim: %d\n",IC.altitude,max,min,bad_trim);
+ printf("Alpha: %g, Throttle_pct: %g, Long_control: %g\n\n",Alpha*RAD_TO_DEG,Throttle_pct,Long_control);
+
+ bad_trim=trim_long(200,IC);
+
+ if(bad_trim == 1 )
+ max=IC.altitude;
+ else
+ min=IC.altitude;
+ i++;
+ }
+}
+
void find_trim_stall(int kmax,FILE *out,InitialConditions IC)
{
return(1);
}
-
+int wave_stats(float *var,float *var_rate,int N,FILE *out)
+{
+ int Nc,i,Nmaxima;
+ float varmax,slope,intercept,time,ld,zeta,omegad,omegan;
+ float varmaxima[100],vm_times[100];
+ /*adjust N so that any constant slope region at the end is cut off */
+ i=N;
+ while((fabs(var_rate[N]-var_rate[i]) < 0.1) && (i >= 0))
+ {
+ i--;
+ }
+ Nc=N-i;
+ slope=(var[N]-var[Nc])/(N*0.01 - Nc*0.01);
+ intercept=var[N]-slope*N*0.01;
+ printf("\tRotating constant decay out of data using:\n");
+ printf("\tslope: %g, intercept: %g\n",slope,intercept);
+ printf("\tUsing first %d points for dynamic response analysis\n",Nc);
+ varmax=0;
+ Nmaxima=0;i=0;
+ while((i <= Nc) && (i <= 801))
+ {
+
+ fprintf(out,"%g\t%g",i*0.01,var[i]);
+ var[i]-=slope*i*0.01+intercept;
+ /* printf("%g\n",var[i]); */
+ fprintf(out,"\t%g\n",var[i]);
+ if(var[i] > varmax)
+ {
+ varmax=var[i];
+ time=i*0.01;
+
+ }
+ if((var[i-1]*var[i] < 0) && (var[i] > 0))
+ {
+ varmaxima[Nmaxima]=varmax;
+ vm_times[Nmaxima]=time;
+ printf("\t%6.2f: %8.4f\n",vm_times[Nmaxima],varmaxima[Nmaxima]);
+ varmax=0;Nmaxima++;
+
+ }
+
+ i++;
+ }
+ varmaxima[Nmaxima]=varmax;
+ vm_times[Nmaxima]=time;
+ Nmaxima++;
+ if(Nmaxima > 2)
+ {
+ ld=log(varmaxima[1]/varmaxima[2]); //logarithmic decrement
+ zeta=ld/sqrt(4*LS_PI*LS_PI +ld*ld); //damping ratio
+ omegad=1/(vm_times[2]-vm_times[1]); //damped natural frequency Hz
+ if(zeta < 1)
+ {
+ omegan=omegad/sqrt(1-zeta*zeta); //natural frequency Hz
+ }
+ printf("\tDamping Ratio: %g\n",zeta);
+ printf("\tDamped Freqency: %g Hz\n\tNatural Freqency: %g Hz\n",omegad,omegan);
+ }
+ else
+ printf("\tNot enough points to take log decrement\n");
+/* printf("w: %g, u: %g, q: %g\n",W_body,U_body,Q_body);
+ */
+ return 1;
+}
// Run an iteration of the EOM (equations of motion)
int main(int argc, char *argv[]) {
double save_alt = 0.0;
- int multiloop=1,k=0,i,j;
+ int multiloop=1,k=0,i,j,touchdown,N;
double time=0,elev_trim,elev_trim_save,elevator,speed,cmcl;
FILE *out;
double hgain,hdiffgain,herr,herrprev,herr_diff,htarget;
SCALAR *control[7];
SCALAR *state[7];
float old_state,effectiveness,tol,delta_state,lctrim;
+ float newcm,lastcm,cmalpha,td_vspeed,td_time,stop_time;
+ float h[801],hdot[801],altmin,lastAlt,theta[800],theta_dot[800];
if(argc < 6)
{
IC.latitude=47.5299892; //BFI
IC.longitude=122.3019561;
Runway_altitude = 18.0;
+
IC.altitude=strtod(argv[2],NULL);
+ printf("h: %g, argv[2]: %s\n",IC.altitude,argv[2]);
IC.vc=strtod(argv[1],NULL);
IC.alpha=0;
IC.beta=0;
- IC.gamma=strtod(argv[3],NULL);
- IC.use_gamma_tmg=1;
+ IC.theta=strtod(argv[3],NULL);
+ IC.use_gamma_tmg=0;
IC.phi=0;
- IC.psi=10;
- IC.weight=2300;
+ IC.psi=0;
+ IC.weight=2400;
IC.cg=0.25;
- IC.flap_handle=0;
- IC.long_control=strtod(argv[4],NULL);
+ IC.flap_handle=10;
+ IC.long_control=0;
IC.rudder_pedal=0;
- printf("IC.vc: %g\n",IC.vc);
+
ls_ForceAltitude(IC.altitude);
fgLaRCsimInit(0.01);
- printf("\nLong_control: %g\n\n",Long_control);
+ setIC(IC);
+ printf("Dx_cg: %g\n",Dx_cg);
+ V_down=strtod(argv[4],NULL);;
+ ls_loop(0,-1);
+ i=0;time=0;
+ IC.long_control=0;
+ altmin=Altitude;
+ printf("\tAltitude: %g, Theta: %g, V_down: %g\n\n",Altitude,Theta*RAD_TO_DEG,V_down);
+
+ while(time < 5.0)
+ {
+ printf("Time: %g, Flap_handle: %g, Flap_position: %g, Transit: %d\n",time,Flap_handle,Flap_Position,Flaps_In_Transit);
+ if(time > 2.5)
+ Flap_handle=20;
+ else if (time > 0.5)
+ Flap_handle=20;
+ ls_update(1);
+ time+=0.01;
+ }
- IC.altitude=1000;
- setIC(IC);
- ls_loop(0.0,-1);
- IC.flap_handle=10;
- setIC(IC);
- ls_loop(0.0,-1);
- IC.flap_handle=20;
- setIC(IC);
- ls_loop(0.0,-1);
- IC.flap_handle=30;
- setIC(IC);
- ls_loop(0.0,-1);
- /* find_trim_stall(200,out,IC);
+ /*out=fopen("drop.out","w");
+ N=800;touchdown=0;
- IC.vc=120;
- IC.altitude=8000;
- IC.weight=2300;
- IC.cg=0.25;
- IC.flap_handle=0;
-
+ while(i <= N)
+ {
+ ls_update(1);
+ printf("\tAltitude: %g, Theta: %g, V_down: %g\n\n",D_cg_above_rwy,Theta*RAD_TO_DEG,V_down);
+ fprintf(out,"%g\t%g\t%g\t%g\t%g\t%g\n",time,D_cg_above_rwy,Theta*RAD_TO_DEG,V_down,F_Z_gear/1000.0,V_rel_ground);
+ h[i]=D_cg_above_rwy;hdot[i]=V_down;
+ theta[i]=Theta; theta_dot[i]=Theta_dot;
+ if(D_cg_above_rwy < altmin)
+ altmin=D_cg_above_rwy;
+ if((F_Z_gear < -10) && (! touchdown))
+ {
+ touchdown=1;
+ td_vspeed=V_down;
+ td_time=time;
+ }
+ time+=0.01;
+ i++;
+ }
+ while(V_rel_ground > 1)
+ {
+ if(Brake_pct < 1)
+ {
+ Brake_pct+=0.02;
+ }
+ ls_update(1);
+ time=i*0.01;
+ fprintf(out,"%g\t%g\t%g\t%g\t%g\t%g\n",time,D_cg_above_rwy,Theta*RAD_TO_DEG,V_down,F_Z_gear/1000.0,V_rel_ground);
+ i++;
+ }
+ stop_time=time;
+ while((time-stop_time) < 5.0)
+ {
+ ls_update(1);
+ time=i*0.01;
+ fprintf(out,"%g\t%g\t%g\t%g\t%g\t%g\n",time,D_cg_above_rwy,Theta*RAD_TO_DEG,V_down,F_Z_gear/1000.0,V_rel_ground);
+ i++;
+ }
+ fclose(out);
- setIC(IC);
- printIC(IC);
- k=trim_long(100,IC);
-
- printf("Flap_handle: %g, Flap_Position: %g\n",Flap_handle,Flap_Position);
- printf("k: %d, %g knots, %g lbs, %g %%MAC\n",k,V_calibrated_kts,Weight,Cg);
- printf("wdot: %g, udot: %g, qdot: %g\n",W_dot_body,U_dot_body,Q_dot_body);
- printf("Alpha: %g, Throttle_pct: %g, Long_control: %g\n\n",Alpha,Throttle_pct,Long_control);
-
- printf("Cme: %g, elevator: %g, Cmde: %g\n",elevator*Cmde,elevator,Cmde);
-
- IC.cg=0.155;
- setIC(IC);
- k=trim_long(100,IC);
+ printf("Min Altitude: %g, Final Alitutde: %g, Delta: %g\n",altmin, h[N], D_cg_above_rwy-altmin);
+ printf("Vertical Speed at touchdown: %g, Time at touchdown: %g\n",td_vspeed,td_time);
+ printf("\nAltitude response:\n");
+ out=fopen("alt.out","w");
+ wave_stats(h,hdot,N,out);
+ fclose(out);
+ out=fopen("theta.out","w");
+ printf("\nPitch Attitude response:\n");
+ wave_stats(theta,theta_dot,N,out);
+ fclose(out);*/
- printf("Flap_handle: %g, Flap_Position: %g\n",Flap_handle,Flap_Position);
- printf("k: %d, %g knots, %g lbs, %g %%MAC\n",k,V_calibrated_kts,Weight,Cg);
- printf("wdot: %g, udot: %g, qdot: %g\n",W_dot_body,U_dot_body,Q_dot_body);
- printf("Alpha: %g, Throttle_pct: %g, Long_control: %g\n\n",Alpha,Throttle_pct,Long_control);
- printf("Cme: %g, elevator: %g, Cmde: %g\n",elevator*Cmde,elevator,Cmde);
-
- IC.cg=0.364;
- setIC(IC);
- k=trim_long(100,IC);
- printf("Flap_handle: %g, Flap_Position: %g\n",Flap_handle,Flap_Position);
+ /*printf("Flap_handle: %g, Flap_Position: %g\n",Flap_handle,Flap_Position);
printf("k: %d, %g knots, %g lbs, %g %%MAC\n",k,V_calibrated_kts,Weight,Cg);
printf("wdot: %g, udot: %g, qdot: %g\n",W_dot_body,U_dot_body,Q_dot_body);
printf("Alpha: %g, Throttle_pct: %g, Long_control: %g\n\n",Alpha,Throttle_pct,Long_control);
- /* do_trims(400,out,IC); */
+
/* ls_loop(0.0,-1);
projects / flightgear.git / blobdiff