1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
|
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* Déclaration des variables du problèmes */
const double m=3;
const double p=1.22;
const double S=0.3;
const double alpha0=10;
const double k=980;
const double Pi=3.141592654;
int vent; /*Permet de choisir le type de vent*/
double coeffw; /*coefficient devant w*/
double alpha;
/* Convertit un angle en degré */
double rad2deg(double angle)
{
return angle*180/Pi;
}
/* Convertit un angle en radian */
double deg2rad(double angle)
{
return angle*Pi/180;
}
/* fonction créant un vent sinusoidal */
double ventSinus(double t)
{
double u,w0=sqrt(k/m);
if(t <= 0) u=0;
else u=7.5+1.5*sin(coeffw*w0*t);
return u;
}
/* fonction créant un vent de type échelon */
double ventEchelon(double t)
{
double u;
if(t<1) u=0;
else
{
if((1<=t) && (t<4)) u=1.5;
else u=0;
}
return u;
}
/* calcul en coefficient de portance en fonction de l'angle */
double coeffPortance(double angle)
{
double Cz,a=deg2rad(angle);
if(-18<=angle && angle<=18) Cz=2*Pi*a;
else Cz=0;
return Cz;
}
/* fonction de la 1er eq différentielle */
double f(double z, double v, double t)
{
return v;
}
/* fonction de la 2er eq différentielle */
double g(double z, double v, double t)
{
double Cz, u;
if(vent==1) u=ventSinus(t);
else u=ventEchelon(t);
alpha=alpha0-atan(v/u);
Cz=coeffPortance(alpha);
return Cz*p*S*u*sqrt(u*u+v*v)/(2*m)-k*z/m;
}
typedef double *fonc(double z, double v, double t);
/* fonction de Runge-Kutta utilisant les fonctions f et g */
/* et permettent de définir z(n+1) en fonction de z(n) */
/* idem pour u(n+1) */
void rk4(double *z, double *v, double *t, double dt, fonc f, fonc g)
{
double kz1=*f(*z,*v,*t)*dt;
double kv1=*g(*z,*v,*t)*dt;
double kz2=*f(*z+kz1/2,*v+kv1/2,*t+dt/2)*dt;
double kv2=*g(*z+kz1/2,*v+kv1/2,*t+dt/2)*dt;
double kz3=*f(*z+kz2/2,*v+kv2/2,*t+dt/2)*dt;
double kv3=*g(*z+kz2/2,*v+kv2/2,*t+dt/2)*dt;
double kz4=*f(*z+kz3,*v+kv3,*t+dt)*dt;
double kv4=*g(*z+kz3,*v+kv3,*t+dt)*dt;
*z=*z+(kz1+2*kz2+2*kz3+kz4)/6;
*v=*v+(kv1+2*kv2+2*kv3+kv4)/6;
}
/* programme principal */
int main()
{
FILE *fp;
double *pz,*pv,*pt;
double z,v,t,dt=0.1;
fonc pf,pg;
z=0;v=0;
pz=&z;pv=&v;pt=&t;
printf("type de vent :\n 1 : vent sinusoidal\n 2 : vent indiciel\nChoix ? ");
scanf("%d",&vent);
if(vent!=1 && vent!=2){
printf("erreur");
return 0;
}
printf("\nValeur de w :\nw=k.w0 avec k=");
scanf("%lf",&coeffw);
/* Appel de rk4 avec les bons paramètres.... */
system("PAUSE");
return 0;
} |
Partager