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;
} |
