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 128 129 130 131 132 133
| #include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define PI 3.14159265359
#define R 8.144621
#define T 298.15
#define F 96485
Q
#define E_IMP 0.6
#define ALPHA 0.5
#define K1 0.1
#define K2 2E-7
#define TAU 1.0
#define N_MAX 60
#define S 1.96E-7
#define C 0.01
#define D 0.00001
/// calcul de R(x)
double r(int x)
{
return pow(x+1, (3.0/2)) - 2*pow(x, (3.0/2)) + pow(x-1, (3.0/2));
}
/// calcul de psie etoile (n) avec mémoire des valeurs déjà calculées
double psie(int n)
{
static const double dt = 1.0 / N_MAX;
static const double ksi = -ALPHA*F/(R*T) * (E_IMP-K1-0.5*log(TAU));
static const double common = 4.0/3*pow(PI, (-1.0/2))*pow(1.0/N_MAX, (1.0/2));
if(n==1)
{
return pow(PI*dt,(-1.0/2))*1.0/(1+exp(-ksi));
}
else
{
double sum = 0;
int j;
for (j=1;j<n;++j)
{
sum += psie(j)*r(n-j);
}
return(1+common*sum)/(exp(-ksi)+common);
}
static double mem[N_MAX];
mem[0] = pow(PI*dt,(-1.0/2))*1.0/(1+exp(-ksi));
if(mem[n-1] == 0)
{
double sum = 0;
int j;
for (j=1;j<n;++j)
{
sum += psie(j)*r(n-j);
}
mem[n-1] = (1+common*sum)/(exp(-ksi)+common);
}
return mem[n-1];
}
/// calcul de omega(n) avec mémoire des valeurs déjà calculées
double omega(int n)
{
static const double dt = 1.0 / N_MAX;
static const double p = K2 * pow(TAU,1.0/2);
static double mem[N_MAX];
mem[0] = (p*dt*psie(n=1))/(2+p*dt*psie(n=1));
if(mem[n-1] == 0)
{
double sum = 0;
int j;
for (j=1;j<n;++j)
{
sum += (1-omega(j))*psie(j);
}
mem[n-1] = (p*dt*psie(n)+sum)/(2+p *dt*psie(n));
}
return mem[n-1];
}
/// calcul de psi(n) avec mémoire des valeurs déjà calculées
double psi(int n)
{
return (1-omega(n=N_MAX))*psie(n=N_MAX);
}
/// calcul de i avec mémoire des valeurs déjà calculées
double i(int n)
{
double sum = 0;
double t;
static double mem[N_MAX];
mem[0] = 0;
{
mem[n-1] = psi(n=N_MAX)*F*S*C*pow(D,1.0/2)/pow((TAU),1.0/2);
}
return mem[n-1];
}
int main()
{
/// création du nom de fichier avec tous les paramètres
//char* filename = "resultat.csv";
char *filename = malloc(255*sizeof(char));
sprintf(filename,"filename_%d_%1.1f_%1.1f_%1.1f_%1.1f", N_MAX, ALPHA, E_IMP, TAU, K1);
strcat(filename, ".csv");
/// verifier que le fichier est bien ouvert
FILE *fp = fopen(filename, "w");
if(fp == NULL) {
printf("Cannot open file.\n");
exit(1);
}
/// on fait tous les calculs
int j;
for (j=1; j<=N_MAX; ++j)
{
//printf("psi(%d) = %e\n", i, psi(i));
//fprintf(fp, "%e\n", psi(i));
printf("psie,omega(%d|%f),i = %1.2e,%1.2e,%1.2e\n", (j*TAU/N_MAX), psie(j), omega(j),i);
fprintf(fp, "%1.2e,%1.2e,%1.2e\n", (j*TAU/N_MAX), psie(j), omega(j),i);
}
fclose(fp);
return 0;
} |
Partager