TheBoussinesqModel  3.2.1
 All Data Structures Files Functions Variables Typedefs Macros Pages
geo_statistic.09375.c
Go to the documentation of this file.
1 
2 /* MATH2 CONTAINS ALGEBRAIC ROUTINES FOR GEOtop AND OTHER MODELS
3 MATH2 Version 0.9375 KMackenzie
4 
5 file geo_statistic.09375.c
6 
7 Copyright, 2009 Stefano Endrizzi, Emanuele Cordano, Matteo Dall'Amico and Riccardo Rigon
8 
9 This file is part of MATH2.
10  MATH2 is free software: you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  MATH2 is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with this program. If not, see <http://www.gnu.org/licenses/>.
22 */
23 
24 
25 
26 #include "turtle.h"
27 #include "t_datamanipulation.h"
28 #include "geo_statistic.09375.h"
29 #include "linearalgebra.h"
30 
31 /*void ordi_kriging(DOUBLEMATRIX *pesi,DOUBLEMATRIX *coord,DOUBLEMATRIX *Z0,DOUBLEVECTOR *U,
32  DOUBLEVECTOR *V, double *scala_integr1, double *varianza1);
33 
34 void variogramma(DOUBLEMATRIX *coord,DOUBLEMATRIX *Z0,DOUBLEMATRIX *variogr,DOUBLEVECTOR *U,
35  double scala_integr, double varianza);
36 
37 double gamma1(double r, double scala_integr, double varianza);
38 
39 int ludcmp(SHORTVECTOR *indx, DOUBLEMATRIX *var);
40 int lubksb(DOUBLEMATRIX *var, SHORTVECTOR *indx,DOUBLEVECTOR *gam);
41 */
42 
78 void ordi_kriging2(DOUBLEMATRIX *weights, DOUBLEVECTOR *E, DOUBLEVECTOR *N, DOUBLEMATRIX *Z0, T_INIT *UV, double int_scale, double variance)
79 
80 {
81 
82 long meno,i,j,k,row,col,n,ii,ncol;
83 short mode;
84 double r,x,y,dx,dy,coordx,coordy;
85 DOUBLEMATRIX *var,*variogr;
86 DOUBLEVECTOR *gam;
87 SHORTVECTOR *indx;
88 
89 /* coordinate vertice basso sx matrice */
90 /* controllare la convenzione FluidTurtle!
91 in r.out.turtle mi sembra essere: {dx,dy,S,W}
92 */
93 
94 coordy=UV->U->co[3];
95 coordx=UV->U->co[4];
96 
97 mode=0;
98 meno=0;
99 row=Z0->nrh;
100 col=Z0->nch;
101 ncol=row*col;
102 dx=UV->U->co[1];
103 dy=UV->U->co[2];
104 n=E->nh; /* numero stazioni */
105 
106 
107 if(n>1){ /* nn>1: piu stazioni: faccio il kriging */
108  var=new_doublematrix(n+1,n+1);
109  gam=new_doublevector(n+1);
110  indx=new_shortvector(n+1);
111  for(i=1;i<=n+1;i++){
112  indx->co[i]=0;
113  }
114  /* Applico la funzione del kriging*/
115  /* valuto la coordinata di ogni stazione rispetto alla coordinata dell'angolo
116  sinistro in basso della mia matrice*/
117  /* Calcolo la funzione var[i][j] definita dalla distanza tra la stazione[i]
118  rispetto alla stazione j*/
119  for(i=1;i<=n;i++){
120  for(j=1;j<=n;j++){
121  x=E->co[i]-E->co[j];
122  y=N->co[i]-N->co[j];
123  r=sqrt(pow(x,(double)2)+pow(y,(double)2)); /* r: distanza */
124  /************* calcolo il variogramma con gamma1(r,int_scale,variance) in geo_statistic.c *********/
125  var->co[i][j]=gamma2(r, int_scale, variance);
126  /**************************************************************/
127  }
128  }
129  for(i=1;i<=n;i++){
130  var->co[i][n+1]=1;
131  }
132  for(j=1;j<=n;j++){
133  var->co[n+1][j]=1;
134  }
135  var->co[n+1][n+1]=0;
136  /*Calcolo le coordinate di ogni pixel rispetto allo stesso punto di riferimento*/
137  /* variogr: doublematrix (nch*nrh, n stazioni) */
138  variogr=new_doublematrix(ncol,n);
139  /************* chiamo variogramma in geo_statistic.c *********/
140  /* calcola la variogramma tra i pixel e le stazioni */
141  variogramma2(E, N , Z0, variogr, UV->U, int_scale, variance);
142  /**************************************************************/
143  if(n>1){
144  for(i=1;i<=row;i++){
145  for(j=1;j<=col;j++){
146  /*Calcolo la distanza di ogni pixel rispetto a tutte le stazioni*/
147  ii=(i-1)*col+j;
148  for(k=1;k<=n;k++){
149  /* gam: vettore con la distanza del pixel dalle stazioni */
150  gam->co[k]=variogr->co[ii][k];
151  }
152  gam->co[n+1]=1;
153  /*Le seguenti subroutine mi risolvono il sistema che avra' come matrice gli elementi
154  var relativi alle stazioni, come vettore dei termini noti il vettore gam relativo
155  ad ogni pixel. La prima routine mi riporta la matrice var modificata per eseguire una
156  sotituzione indietro, questo e' sufficiente eseguirla un unica volta */
157  if(i==1 && j==1){
158  /************* chiamo ludcmp in linearalgebra.c *********/
159  ludcmp(indx,var);
160  /**************************************************************/
161  }
162  /************* chiamo lubksb in linearalgebra.c *********/
163  lubksb(var,indx,gam);
164  /**************************************************************/
165  /* metto i risultati ottenuti in gam nella matrice pesi */
166  if(Z0->co[i][j]!=UV->V->co[2]){
167  for(k=1;k<=n;k++){
168  weights->co[ii][k]=gam->co[k];
169  }
170  }
171  }
172  }
173  }else{
174  for(i=1;i<=row;i++){
175  for(j=1;j<=col;j++){
176  /*Calcolo la distanza di ogni pixel rispetto a tutte le stazioni*/
177  ii=(i-1)*col+j;
178  weights->co[ii][1]=variogr->co[ii][2]/(variogr->co[ii][1]+variogr->co[ii][2]);
179  weights->co[ii][2]=variogr->co[ii][1]/(variogr->co[ii][1]+variogr->co[ii][2]);
180  }
181  }
182  }
183  free_doublevector(gam);
184  free_doublematrix(var);
185  free_doublematrix(variogr);
186  free_shortvector(indx);
187 }else{ /* n=1: 1 stazione: non faccio il kriging e pongo tutti i pesi=1 */
188  //printf("THERE IS ONE STATION\n");
189  for(i=1;i<=Z0->nrh;i++){
190  for(j=1;j<=Z0->nch;j++){
191  ii=(i-1)*col+j;
192  if(Z0->co[i][j]!=UV->V->co[2]){
193  weights->co[ii][1]=1;
194  }
195  }
196  }
197 }
198 /* tolgo eventuali pesi < 0 che non hanno senso */
199 for(i=1;i<=ncol;i++){
200  for(j=1;j<=n;j++){
201  if(weights->co[i][j]<=0) weights->co[i][j]=0;
202  }
203 }
204 }
205 
206 
243 void variogramma2(DOUBLEVECTOR *E, DOUBLEVECTOR *N, DOUBLEMATRIX *Z0, DOUBLEMATRIX *variogr, DOUBLEVECTOR *U, double int_scale, double variance)
244 
245 {
246 long ii,i,j,k;
247 double X,Y,rx,ry,rz,h2;
248 
249 ii=0;
250 for(i=1;i<=Z0->nrh;i++){
251  Y=U->co[3]+(Z0->nrh-i)*U->co[2]+(U->co[2]/(double)2);
252  for(j=1;j<=Z0->nch;j++){
253  ii++;
254  for(k=1;k<=E->nh;k++){
255  X=U->co[4]+(j-1)*U->co[1]+(U->co[1]/(double)2);
256  rx=X-E->co[k];
257  ry=Y-N->co[k];
258  rz=0;
259  h2=sqrt(pow(rx,(double)2)+pow(ry,(double)2)+pow(rz,(double)2));
260  if(h2!=0)variogr->co[ii][k]=gamma2(h2,int_scale,variance);
261  }
262  }
263 }
264 }
265 
266 
294 double gamma2(double r, double scala_integr, double varianza)
295 
296 {
297 double gamma2;
298 
299 gamma2=varianza*(1.0-exp(-r/scala_integr));
300 
301 return gamma2;
302 }
303