computeIjxy.c 3.25 KB
Newer Older
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
#include "stdio.h"
#include "string.h"
#include <stdlib.h>
#include <sys/dir.h>
#include <sys/stat.h>
#include <sys/types.h>
#include "math.h"

#define NB_PT_SAFRAN 9892
#define NB_TIME_SAFRAN 8760

#define FIRST_YEARS 2011
#define LAST_YEARS 2012
#define T_ZERO_OBSOLU -273.15
#define PATH_TO_DATA_SAFRAN "/home/olivierb/solvers/trunk/SandBox/FF/AedesAlbopictus/intraAnnuel/donnees_SAFRAN/AUTOROUTES/DATA/"

//11.128007398517115   0.078865315431034  12.651026892598505 42.345187120388388  -0.145204446756075

double psi=11.128007398517115;
double alpha=0.078865315431034;
double delta=12.651026892598505;
double Tmax=42.345187120388388;
double lambda=-0.145204446756075 ;


//interval de temps:
/*
on veut les dates comprise entre le 1 mars et le 31 dec
Les donnees de 2011 commence le 1aout 2011.
on garde les dates du 1aout au 31 dec : de 0 a 3666
on ignore les dates du 1j au 1 mars ie de 3666 a 59*24+3666=5082
 */
double Tmin=0;
double Vmin;
double Vmax=-100;
double rhofct(double T){
  double res=psi*(exp(alpha*T)-exp(alpha*Tmax-(Tmax-T)/delta))+lambda;
  if (T<Tmin){
    Tmin=T;
    Vmin=res;
  }
  if (Vmax<res)
    Vmax=res;
  return res;
}

#define PERIOD_HEURE 24

void safranToIxyTday( ){
  char inputDataSafran[512];
  double *pValT=(double*)malloc(NB_PT_SAFRAN*sizeof(double));
  
  int step=1;
  double * pTday=(double *)calloc(NB_PT_SAFRAN,sizeof(double));
  double * pIjxy=(double *)calloc(NB_PT_SAFRAN,sizeof(double));
  FILE * foutput=fopen("IjxyTj.txt","w");
  int cmp=0;
  for (int numA=FIRST_YEARS;numA<LAST_YEARS;numA++){
    int k;
    for (int step=1; step<NB_TIME_SAFRAN+1; step++){
       if (step > 3666 && step < 5082){
         cmp=0;
         continue;
       }
      //lecture des donnees SAFRAN, ie raster incomplet 
      //sprintf(filename,"%s%i",pathData,i);
      sprintf(inputDataSafran,"%s/%i/TAIR_%i",PATH_TO_DATA_SAFRAN,numA,step);
      FILE *fTempSafran=fopen(inputDataSafran,"rb");
      if (!fTempSafran)
        continue;
      int rfread=fread(pValT, sizeof(double), NB_PT_SAFRAN, fTempSafran);
      if (rfread!=NB_PT_SAFRAN){
        printf("can't read %i double in file %s, read only %i\n",NB_PT_SAFRAN,inputDataSafran,rfread);
        return;
      }
      fclose(fTempSafran);
      cmp++;
      for(int locSafran=0;locSafran<NB_PT_SAFRAN;locSafran++){
        pValT[locSafran]+=T_ZERO_OBSOLU;
        pTday[locSafran]+=pValT[locSafran];
        pIjxy[locSafran]+=rhofct(pValT[locSafran]);
      }
      if (!(step%PERIOD_HEURE)){
          for(int locSafran=0;locSafran<NB_PT_SAFRAN;locSafran++){
            //if (pTday[locSafran]/24>10)
            //  printf("%e %e\n",pTday[locSafran]/24,pIjxy[locSafran]);
            if (cmp==PERIOD_HEURE){
              fprintf(foutput,"%lf %lf\n",pTday[locSafran]/(1.0*PERIOD_HEURE),pIjxy[locSafran]/(1.0*PERIOD_HEURE));
              if (pIjxy[locSafran]/(1.0*PERIOD_HEURE)>Vmax)
                printf("BUG %i\n",step);
            }
            pTday[locSafran]=0;
            pIjxy[locSafran]=0;
          }
          cmp=0;
      }
    }
  }
  free(pTday);
  free(pIjxy);
  free(pValT);
  fclose(foutput);
}

int main(){
  FILE *pf=fopen("rho.txt","w");
  for(int ii=-20;ii<45;ii++)
    fprintf(pf,"%e %e\n",1.0*ii,rhofct(1.0*ii));
  fclose(pf);
  
  
  safranToIxyTday();
  printf("%e %e %e\n",Tmin,Vmin,Vmax);
  return 0;
}