source: issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp@ 18446

Last change on this file since 18446 was 18446, checked in by abuzzi, 11 years ago

MODIF: final version not tested
File size: 19.9 KB
RevLine 
[18434]1#include <iostream>
2#include <cmath>
3#include <fstream>
4#include <string>
5#include <cstdio>
6#include <cstdlib>
7#include <time.h>
8#include <cassert>
9#include <gsl/gsl_multifit.h>
10
11using namespace std;
12
13class Matrix{/*{{{*/
14 private:
15 int M; /*Number of lines */
16 int N; /*Number if Columns */
17 double *values;
18 public:
19 Matrix(int m_in,int n_in){/*{{{*/
20 this->M = m_in;
21 this->N = n_in;
22 this->values = new double[M*N];
23 }/*}}}*/
24 ~Matrix(){/*{{{*/
25 delete [] this->values;
26 }/*}}}*/
27 void Echo(void){/*{{{*/
28 for(int i=0;i<M;i++){
29 for(int j=0;j<N;j++){
30 cout << " " << this->values[i*N+j];
31 }
32 cout << endl;
33 }
34 }/*}}}*/
35 void SetValue(int i,int j,double value){/*{{{*/
36 this->values[i*N+j] = value;
37 }/*}}}*/
38 double GetValue(int i,int j){/*{{{*/
39 return this->values[i*N+j];
40 }/*}}}*/
41 void GetSize(int* pM,int* pN){/*{{{*/
42 *pM = this->M;
43 *pN = this->N;
44 }/*}}}*/
45 double* GetPointer(void){/*{{{*/
46 return this->values;
47 }/*}}}*/
48 void MatrixSum(Matrix* A,Matrix* B){/*{{{*/
49 /*Check that sizes are compatible*/
50 int M_B,N_B,M_A,N_A;
51 B->GetSize(&M_B,&N_B);
52 A->GetSize(&M_A,&N_A);
53 assert(this->M==M_B && this->N==N_B);
54 assert(this->M==M_A && this->N==N_A);
55 for(int i=0;i<this->M;i++){
56 for(int j=0;j<this->N;j++){
57 this->SetValue(i,j,A->GetValue(i,j) + B->GetValue(i,j));
58 }
59 }
60 }/*}}}*/
61 void MatrixDiff(Matrix* A,Matrix* B){/*{{{*/
62 /*Check that sizes are compatible*/
63 int M_B,N_B,M_A,N_A;
64 B->GetSize(&M_B,&N_B);
65 A->GetSize(&M_A,&N_A);
66 assert(this->M==M_B && this->N==N_B);
67 assert(this->M==M_A && this->N==N_A);
68 for(int i=0;i<this->M;i++){
69 for(int j=0;j<this->N;j++){
70 this->SetValue(i,j,A->GetValue(i,j) - B->GetValue(i,j));
71 }
72 }
73 }/*}}}*/
74 void MatrixAbs(Matrix* A){/*{{{*/
75 for(int i=0;i<this->M;i++){
76 for(int j=0;j<this->N;j++){
77 this->SetValue(i,j,fabs(A->GetValue(i,j)));
78 }
79 }
80 }/*}}}*/
81 void MatrixEqual(Matrix* A){/*{{{*/
82 for(int i=0;i<this->M;i++){
83 for(int j=0;j<this->N;j++){
84 this->SetValue(i,j,A->GetValue(i,j));
85 }
86 }
87 }/*}}}*/
88 double MatrixInternSum(){/*{{{*/
89 double sum=0;
90 for(int i=0;i<this->M;i++){
91 for(int j=0;j<this->N;j++){
92 sum+=this->GetValue(i,j);
93 }
94 }
95 return sum;
96 }/*}}}*/
97 void ExtractLine(Matrix* A,int i){/*{{{*/
98 /* Check that the size of A is compatible */
99 int M_A,N_A;
100 A->GetSize(&M_A,&N_A);
101 assert(M_A==1 && this->N==N_A);
102 for(int j=0;j<this->N;j++){
103 A->SetValue(0,j,this->GetValue(i,j));
104 }
105 }/*}}}*/
106 void ExtractColumn(Matrix* A,int j){/*{{{*/
107 /* Check that the size of A is compatible */
108 int M_A,N_A;
109 A->GetSize(&M_A,&N_A);
110 assert(N_A==1 && this->M==M_A);
111 for(int i=0;i<this->M;i++){
112 A->SetValue(i,0,this->GetValue(i,j));
113 }
114 }/*}}}*/
115 void AddNumber(double a){/*{{{*/
116 for(int i=0;i<this->M;i++){
117 for(int j=0;j<this->N;j++){
118 this->SetValue(i,j,this->GetValue(i,j) + a);
119 }
120 }
121 }/*}}}*/
122};/*}}}*/
123
124/*Local prototypes{{{*/
125void makep(Matrix *Pobs,int nx,int ny, int dx, int dy);
126void vec2grid(Matrix *V,Matrix *V1,Matrix *V2,int nx, int ny);
127void msplit( Matrix *m, Matrix *m1,Matrix *m2,double dlevel);
128void plouff(Matrix *g,Matrix *Pobs,Matrix *Pp,Matrix * mesh,Matrix *rho,int dx,int dy, int dn,int m,int n,int l);
129void vec2gridsimple(Matrix *V,Matrix *V1,int nx, int ny);
130void reshape(Matrix* V,Matrix* V1,int nx,int ny);
131double misfit(Matrix* m0,Matrix* evalid,Matrix* gobs,int dlevel,Matrix* Pobs,Matrix* xobs,Matrix* yobs,Matrix* Pp,Matrix* rho1, Matrix* rho2,int dx,int dy,int dn,int nx,int ny, int mx,int my);
132void GSLsquarefit(Matrix** pX,Matrix* A,Matrix* B);
[18435]133double signe(double a);
134void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my);
[18439]135void newmodelgen(Matrix* m0,Matrix* m1,Matrix* bathy,Matrix* icethick,int mx,int my,double T,double ptval,double mmax,double mmax2,double ctr,double sd);
[18446]136double coolshed(double T0,int k,double c,double D);
[18434]137/*}}}*/
138
139int main(){/*{{{*/
140
141 /* Seed the random number generator {{{*/
142 srand (time(NULL)); /*}}}*/
143 /* Define the variables {{{*/
144
145 int dx = 500; /* prism dimension in x-direction */
146 int dy = 500; /* prism dimension in y-direction */
147 int mx = 39; /* number of prisms in x-direction */
148 int my = 60; /* number of prisms in y-direction */
149 int nx = 39; /* number of data points in x-direction */
150 int ny = 60; /* number of data points in y-direction */
151 int dn = 15000; /* distance for neighbouting prisms for gravity calculation */
152 double ptval = 100.; /* max. amount to perturb model */
153 double ptval2 = 100.;
154
155 Matrix *Pobs=new Matrix(nx*ny,2); /* data positions */
156 makep(Pobs,nx,ny,dx,dy);
157// Pobs->Echo();
158
159
160 Matrix *Pp=new Matrix(mx*my,2); /* data positions */
161 makep(Pp,mx,my,dx,dy);
162// Pp->Echo();
163
164 double rhoi = 890; /* ice density */
165 double rhow = 1030; /* water density */
166// double rhos = 2013; /* sediment density */
167 double rhoc = 2670; /* bedrock density */
168
169 Matrix *Rho = new Matrix(1,2);
170 Rho->SetValue(0,0,rhoi);
171 Rho->SetValue(0,1,rhow);
172 Matrix *rho1 = new Matrix(1,3);
173 rho1->SetValue(0,0,rhoi);
174 rho1->SetValue(0,1,rhow);
175 rho1->SetValue(0,2,rhoc);
176 Matrix *rho2 = new Matrix(1,2);
177 rho2->SetValue(0,0,rhoi-rhoc);
178 rho2->SetValue(0,1,rhow-rhoc);
179
180 double ctr=1; /* parameter for filtering */
181 double sd=0.1;
182
183 Matrix *xobs= new Matrix(ny,nx);
184 Matrix *yobs= new Matrix(ny,nx);
185
186 vec2grid(Pobs,xobs,yobs,nx,ny);
187// xobs->Echo();
188// yobs->Echo();
189
190 double dlevel= 860; /* mean level of data acquisition */
191
192 double mmax = 1000; /* max value for layer interfaces */
193 double mmax2 = 1000;
194 double mmax3 = 1000;
195
196 /* control parameter for temperature schedule */
197
198 double ca=0.9; /* for acceptance */
199 double cm=0.5; /* for model perturbation */
[18446]200
201 double T0a = 0.1; /* initial temperature for acceptance */
202 double T0m = 0.9; /* initial temperature for model perturbation */
203 double D = 2; /* dimension of the model */
204 int maxconsecrej = 1000; /* max consecutive rejection */
205 int maxsuccess = 100; /* max number of success within one temperature */
206 double T_min = 1e-10; /* stopping temp */
207 double Tred = 1;
208 double E_min = -1000000;
209 double E_exp = 0.0291; /* expected misfit */
210 int maxiter = 10000;
211 int maxtotaliter = 1000000;
212 double Tol = 1e-10; /* tolerance on misfit */
213 int sfreq = 100;
214
[18434]215 /*}}}*/
216 /* load the data {{{*/
217
218 /* Observed gravity anomaly */
219
220 ifstream file1("store_fa500_36s.txt");
221 Matrix * gobs= new Matrix(nx*ny,1);
222 double inputnumber;
223 for(int i=0;i<ny*nx; i++){
224 file1 >> inputnumber;
225 gobs->SetValue(i,0, inputnumber*1e-5);
226 }
227 file1.close();
228// gobs->Echo();
229
230 /* load data about the ice thickness */
231
232 ifstream file2("store_flag_icethick500.txt");
233 Matrix * icethick= new Matrix(mx*my,1);
234 for(int s=0;s<mx*my; s++){
235 file2 >> inputnumber;
236 icethick->SetValue(s,0,inputnumber);
237 }
238 file2.close();
239// icethick->Echo();
240
241 /* load the batimethry data */
242
243 ifstream file3("store_flag_bathy500.txt");
244 Matrix * bathy= new Matrix(mx*my,1);
245 for(int s=0;s<mx*my; s++){
246 file3 >> inputnumber;
247 bathy->SetValue(s,0,inputnumber);
248 }
249 file3.close();
250// bathy->Echo();
251
252 /* id of grid to evaluate misfit */
253
254
255 ifstream file4("store_flag_eval500.txt");
256 Matrix * evalid= new Matrix(mx*my,1);
257 for(int s=0;s<mx*my; s++){
258 file4 >> inputnumber;
259 evalid->SetValue(s,0,inputnumber);
260 }
261 file4.close();
262// evalid->Echo();
263
264 /* initial guess of the model */
265
266 ifstream file5("m0_140114b.txt");
267 Matrix * mesh_ini= new Matrix(mx*my,3);
268 for(int s=0;s<mx*my; s++){
269 for(int j=0;j<3;j++){
270 file5 >> inputnumber;
271 mesh_ini->SetValue(s,j,inputnumber);
272 }
273 }
274 file5.close();
275// mesh_ini->Echo();
276 /*}}}*/
[18446]277 /* VFSA {{{ */
278
279 /* name of the files to save results */
280 std::ofstream savefile1 ("r_140114b.txt");
281 std::ofstream savefile2("m_140114b.txt");
282
283 /* counters initialization */
284 int success = 0;
285 int finished = 0;
286 int consec = 0;
287 double Ta = T0a;
288 double Tm = T0m;
289 int iterT = 0; /* iteration within a T */
290 int total = 0; /* total number of iteration */
291 int totaliter = 0;
292 int msave = 0;
293 double E_new;
294 double E_final;
295 double dE;
296 double P;
297 double rn;
298 Matrix* m_old = new Matrix(mx *my,3);
299 Matrix* m_min = new Matrix(mx *my,3);
300 Matrix* m_new = new Matrix(mx *my,3);
301 m_old->MatrixEqual(mesh_ini);
302 mesh_ini->Echo();
303
304 /* calculate initial misfit */
305 double E_old=misfit(m_old,evalid,gobs,dlevel,Pobs,xobs,yobs,Pp,rho1,rho2,dx,dy,dn,nx,ny,mx,my);
306
307 /* record initial settings */
308 savefile1 << "P "<< "Ta "<< "Tm "<< "Eold "<< "totaliter "<< "Tred "<< endl;
309 savefile1 << "nan "<< Ta<<" "<< Tm<<" "<< E_old<<" "<< totaliter<<" "<< Tred <<" "<< endl;
310 savefile2 << totaliter<< endl;
311 for(int i=0;i<mx*my;i++){
312 savefile2 << m_old->GetValue(i,0)<<" "<< m_old->GetValue(i,1)<<" "<< m_old->GetValue(i,2)<<endl;
313 }
314
315 /* beginning of the loop */
316
317 while(finished==0){
318
319 iterT++;
320 totaliter++;
321
322 /* stop or reduce T */
323 if(iterT>=maxiter || success>maxsuccess){
324 if(Ta<T_min || total>maxtotaliter || fabs(E_old)<=Tol){
325 finished=1;
326 total+=iterT;
327 break;
328 }
329 else{ /* reduce T */
330 Ta=coolshed(T0a,Tred,ca,D);
331 Tm=coolshed(T0m,Tred,cm,D);
332 total+=iterT;
333 iterT=0;
334 success=1;
335 Tred++;
336 consec=0;
337 }
[18434]338 }
[18446]339
340 /* update model and calculate energy */
341
342 newmodelgen(m_old,m_new,bathy,icethick,mx,my,Tm,ptval,mmax,mmax2,ctr,sd); /* new model */
343 E_new=misfit(m_new,evalid,gobs,dlevel,Pobs,xobs,yobs,Pp,rho1,rho2,dx,dy,dn,nx,ny,mx,my); /* new energy */
344 dE=E_new-E_old; /* energy difference */
345
346 /* acceptance probability */
347
348 P=exp(-dE/Ta);
349
350 /* stop if energy is lower than specified minimum */
351 if (E_new<E_min){
352 m_old->MatrixEqual(m_new);
353 E_old=E_new;
354 break;
355 }
356
357 rn=rand()/double (RAND_MAX);
358
359 /* accept new model or not */
360 if(dE<=0){
361 m_old->MatrixEqual(m_new);
362 E_old=E_new;
363 success++;
364 consec=0;
365 savefile1 << P<<" "<< Ta<<" "<< Tm<<" "<< E_old<<" "<< totaliter<<" "<< Tred <<" "<< endl;
366 if(Ta<1e-3){
367 savefile2 << totaliter<< endl;
368 for(int i=0;i<mx*my;i++){
369 savefile2 << m_old->GetValue(i,0)<<" "<< m_old->GetValue(i,1)<<" "<< m_old->GetValue(i,2)<<endl;
370 }
371 }
372 }
373 else{
374 if(P>rn){
375 m_old->MatrixEqual(m_new);
376 E_old=E_new;
377 success++;
378 consec=0;
379 savefile1 << P<<" "<< Ta<<" "<< Tm<<" "<< E_old<<" "<< totaliter<<" "<< Tred <<" "<< endl;
380 if(Ta<1e-3){
381 savefile2 << totaliter<< endl;
382 for(int i=0;i<mx*my;i++){
383 savefile2 << m_old->GetValue(i,0)<<" "<< m_old->GetValue(i,1)<<" "<< m_old->GetValue(i,2)<<endl;
384 }
385 }
386 }
387 else{
388 consec++;
389 }
390 }
[18434]391 }
392
[18446]393 m_min->MatrixEqual(m_old);
394 E_final=E_old;
395 savefile1 << "nan"<<" "<< "nan"<<" "<< "nan"<<" "<< E_final<<" "<< "nan"<<" "<< "nan" <<" "<< endl;
396 savefile2 << " Mesh final"<< endl;
397 for(int i=0;i<mx*my;i++){
398 savefile2 << m_min->GetValue(i,0)<<" "<< m_min->GetValue(i,1)<<" "<< m_min->GetValue(i,2)<<endl;
399 }
400 /*}}}*/
[18434]401 return 0;
402}/*}}}*/
403
404void GSLsquarefit(Matrix** pX,Matrix* A,Matrix* B){/*{{{*/
405
406 /*GSL Matrices and vectors: */
407 int M,N;
408 double chisq;
409
410 /*Get system size*/
411 A->GetSize(&M,&N);
412
413 /*Initialize gsl matrices and vectors: */
414 gsl_matrix* a = gsl_matrix_alloc(M,N);
415 for(int i=0;i<M;i++){
416 for(int j=0;j<N;j++){
417 gsl_matrix_set (a,i,j,A->GetValue(i,j));
418 }
419 }
420 gsl_vector* b = gsl_vector_alloc(M);
421 for(int i=0;i<M;i++){
422 gsl_vector_set(b,i,B->GetValue(i,0));
423 }
424
425 gsl_vector* x = gsl_vector_alloc(N);
426 gsl_matrix* cov = gsl_matrix_alloc(N,N);
427
428 /*Least square fit: */
429 gsl_multifit_linear_workspace* work = gsl_multifit_linear_alloc(M,N);
430 gsl_multifit_linear (a, b, x, cov, &chisq, work);
431 gsl_multifit_linear_free (work);
432
433 /*Clean up and assign output pointer*/
434 Matrix* X = new Matrix(N,1);
435 for(int j=0;j<N;j++){
436 X->SetValue(j,0,gsl_vector_get(x,j));
437 }
438 *pX = X;
439
440 gsl_matrix_free(a);
441 gsl_vector_free(b);
442 gsl_matrix_free(cov);
443
444}/*}}}*/
445void makep(Matrix *Pobs,int nx,int ny, int dx, int dy){/*{{{*/
446 for(int i=0;i<ny;i++){
447 for(int j=0;j<nx;j++){
448 Pobs->SetValue(j+nx*i,0,j*dx);
449 Pobs->SetValue(j+nx*i,1,i*dy);
450 }
451 }
452}/*}}}*/
453void vec2grid(Matrix *V,Matrix *V1,Matrix *V2,int nx, int ny){/*{{{*/
454 for(int i=0;i<ny;i++){
455 for (int j=0;j<nx;j++){
456 V1->SetValue(i,j, V->GetValue(j+nx*i,0));
457 V2->SetValue(i,j, V->GetValue(j+nx*i,1));
458 }
459 }
460}/*}}}*/
461void msplit( Matrix *m, Matrix *m1,Matrix *m2,double dlevel){/*{{{*/
462 int sizem1,sizem2;
463 m->GetSize(&sizem1,&sizem2);
464 for(int i=0;i<sizem1;i++){
465 for(int j=0;j<sizem2+1;j++){
466 if(j<sizem2){
467 m1->SetValue(i,j,1e-10*(sizem2+1-j));
468 m2->SetValue(i,j,m->GetValue(i,j));
469 if(m->GetValue(i,j)<0){
470 m1->SetValue(i,j,m->GetValue(i,j));
471 m2->SetValue(i,j,i*1e-10);
472 }
473 }
474 else{
475 m1->SetValue(i,j,1e-10);
476 }
477 }
478 }
479 m1->AddNumber(dlevel);
480 m2->AddNumber(dlevel);
481}/*}}}*/
482void plouff(Matrix *g,Matrix *Pobs,Matrix *Pp,Matrix * mesh,Matrix *rho,int dx,int dy, int dn,int m,int n,int l){/*{{{*/
483 double gg=6.673e-11;
484 int si,sj,id,s;
485 double R,Q,P;
486 Matrix *xp= new Matrix(1,2);
487 Matrix *yp= new Matrix(1,2);
488 Matrix *xpp= new Matrix(1,2);
489 Matrix *ypp= new Matrix(1,2);
490 Matrix *U= new Matrix(l,4);
491 Matrix *U1=new Matrix(1,4);
492 Matrix *U2=new Matrix(1,4);
493 Matrix *gl= new Matrix(1,l-1);
494 bool test=true;
495 for(int c=0;c<n;c++){
496 g->SetValue(c,0,0);
497 for(int a=0;a<m;a++){
498 test=true;
499 xp->SetValue(0,0,Pp->GetValue(a,0)-Pobs->GetValue(c,0));
500 xp->SetValue(0,1,Pp->GetValue(a,0)-Pobs->GetValue(c,0)+dx);
501 if(xp->GetValue(0,0)<0 && xp->GetValue(0,0)<xp->GetValue(0,1) && xp->GetValue(0,0)*xp->GetValue(0,1)>=0){
502 xpp->SetValue(0,0,xp->GetValue(0,1));
503 xpp->SetValue(0,1,xp->GetValue(0,0));
504 xp->MatrixAbs(xpp);
505 }
506 yp->SetValue(0,0,Pp->GetValue(a,1)-Pobs->GetValue(c,1));
507 yp->SetValue(0,1,Pp->GetValue(a,1)-Pobs->GetValue(c,1)+dy);
508 if(yp->GetValue(0,0)<0 && yp->GetValue(0,0)<yp->GetValue(0,1) && yp->GetValue(0,0)*yp->GetValue(0,1)>=0){
509 ypp->SetValue(0,0,yp->GetValue(0,1));
510 ypp->SetValue(0,1,yp->GetValue(0,0));
511 yp->MatrixAbs(ypp);
512 }
513 P=sqrt(xp->GetValue(0,0)*xp->GetValue(0,0)+yp->GetValue(0,0)*yp->GetValue(0,0));
514 if(P>dn){
515 test=false;
516 for(int i=0;i<l-1;i++){
517 gl->SetValue(0,i,0);
518 }
519 }
520 if(test==true){
521 si=1;
522 sj=1;
523 id=0;
524 for(int i=0;i<2;i++){
525 si*=-1;
526 for(int j=0;j<2;j++){
527 si*=-1;
528 s=si*sj;
529 for(int k=0;k<l;k++){
530 R=sqrt(xp->GetValue(0,i)*xp->GetValue(0,i)+yp->GetValue(0,j)*yp->GetValue(0,j)+mesh->GetValue(a,k)*mesh->GetValue(a,k));
531 Q=atan(xp->GetValue(0,i)*yp->GetValue(0,j)/(mesh->GetValue(a,k)*R));
532 U->SetValue(k,id,s*(mesh->GetValue(a,k)*Q-xp->GetValue(0,i)*log(R+yp->GetValue(0,j))-yp->GetValue(0,j)*log(R+xp->GetValue(0,i))));
533 }
534 id++;
535 }
536 }
537 for(int b=0;b<l-1;b++){
538 U->ExtractLine(U1,b);
539 U->ExtractLine(U2,b+1);
540 gl->SetValue(0,b,rho->GetValue(0,b)*(U1->MatrixInternSum()*(-1)+U2->MatrixInternSum()));
541 }
542 }
543 g->SetValue(c,0,g->GetValue(c,0)+gg*gl->MatrixInternSum());
544 }
545 }
546 delete xp;
547 delete yp;
548 delete xpp;
549 delete ypp;
550 delete gl;
551 delete U;
552 delete U1;
553 delete U2;
554 }/*}}}*/
555void vec2gridsimple(Matrix *V,Matrix *V1,int nx, int ny){/*{{{*/
556 for(int i=0;i<ny;i++){
557 for (int j=0;j<nx;j++){
558 V1->SetValue(i,j, V->GetValue(j+nx*i,0));
559 }
560 }
561}/*}}}*/
562void reshape(Matrix* V,Matrix* V1,int nx,int ny){/*{{{*/
563 for (int i=0;i<ny;i++){
564 for(int j=0;j<nx;j++){
565 V1->SetValue(j+nx*i,0,V->GetValue(i,j));
566 }
567 }
568}/*}}}*/
569double misfit(Matrix* m0,Matrix* evalid,Matrix* gobs,int dlevel,Matrix* Pobs,Matrix* xobs,Matrix* yobs,Matrix* Pp,Matrix* rho1, Matrix* rho2,int dx,int dy,int dn,int nx,int ny, int mx,int my){/*{{{*/
570 Matrix* m1=new Matrix(mx*my,4);
571 Matrix* m2=new Matrix(mx*my,3);
572 Matrix* g1=new Matrix(nx*ny,1);
573 Matrix* g2=new Matrix(nx*ny,1);
574 Matrix* g=new Matrix(nx*ny,1);
575 Matrix* gcalgr=new Matrix(ny,nx);
576 Matrix* gcalvec=new Matrix(mx*my,1);
577 Matrix* df=new Matrix(mx*my,1);
578 Matrix* G=new Matrix(mx*my,3);
579 double a=0;
580 double b=0;
581 double e=0;
582 msplit(m0,m1,m2,dlevel);
583 plouff(g1,Pobs,Pp,m1,rho1,dx,dy,dn,mx*my,nx*ny,4);
584 plouff(g2,Pobs,Pp,m2,rho2,dx,dy,dn,mx*my,nx*ny,3);
585 g->MatrixSum(g1,g2);
586 vec2gridsimple(g,gcalgr,nx,ny);
587 reshape(gcalgr,gcalvec,mx,my);
588 for (int i=0;i<mx*my;i++){
589 df->SetValue(i,0,evalid->GetValue(i,0)*(gobs->GetValue(i,0)-gcalvec->GetValue(i,0)));
590 G->SetValue(i,0,evalid->GetValue(i,0)*Pobs->GetValue(i,0));
591 G->SetValue(i,1,evalid->GetValue(i,0)*Pobs->GetValue(i,1));
592 G->SetValue(i,2,evalid->GetValue(i,0));
593 }
594 Matrix* M = NULL;
595 GSLsquarefit(&M,G,df);
596
597 for (int i=0;i<my;i++){
598 for(int j=0;j<mx;j++){
599 gcalgr->SetValue(i,j,gcalgr->GetValue(i,j)+xobs->GetValue(i,j)*M->GetValue(0,0)+yobs->GetValue(i,j)*M->GetValue(1,0)+M->GetValue(2,0));
600 }
601 }
602 reshape(gcalgr,g,mx,my);
603 for (int i=0;i<mx*my;i++){
604 a=a+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)-g->GetValue(i,0)));
605 b=b+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)+g->GetValue(i,0)));
606 }
607 e=2*a/(a+b);
608 return e;
609}/*}}}*/
[18435]610void newmodelgen(Matrix* m0,Matrix* m1,Matrix* bathy,Matrix* icethick,int mx,int my,double T,double ptval,double mmax,double mmax2,double ctr,double sd){/*{{{*/
611 Matrix* m1gr=new Matrix(my,mx);
612 Matrix* m1grsm=new Matrix(my,mx);
613 Matrix* m1col=new Matrix(mx*my,1);
614 double u=0;
615 double y=0;
616 m1->MatrixEqual(m0);
617 for (int i=0;i<mx*my;i++){
618 if(icethick->GetValue(i,0)==0){
619 u=double (rand())/ double(RAND_MAX);
620 y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
621 m1->SetValue(i,1,m0->GetValue(i,1)+y*ptval);
622 if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
623 m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
624 }
625 if(m1->GetValue(i,1)>=m1->GetValue(i,0)+mmax){
626 m1->SetValue(i,1,m1->GetValue(i,0)+mmax);
627 }
628 }
629 }
630 m1->ExtractColumn(m1col,1);
631 vec2gridsimple(m1col,m1gr,mx,my);
[18439]632 cout<<endl<<endl;
[18435]633 filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
634 reshape(m1grsm,m1col,mx,my);
635 for (int i=0;i<mx*my;i++){
636 if(icethick->GetValue(i,0)==0){
637 m1->SetValue(i,1,m1col->GetValue(i,0));
638 }
639 else{
640 m1->SetValue(i,1,m0->GetValue(i,1));
641 }
642 if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
643 m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
644 }
645 }
646
647 for (int i=0;i<mx*my;i++){
648 if(bathy->GetValue(i,0)==0){
649 u=double (rand())/ double(RAND_MAX);
650 y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
651 m1->SetValue(i,2,m0->GetValue(i,2)+y*ptval);
652 if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
653 m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
654 }
655 if(m1->GetValue(i,2)>=m1->GetValue(i,1)+mmax2){
656 m1->SetValue(i,2,m1->GetValue(i,1)+mmax2);
657 }
658 }
659 }
660 m1->ExtractColumn(m1col,2);
661 vec2gridsimple(m1col,m1gr,mx,my);
662 filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
663 reshape(m1grsm,m1col,mx,my);
664 for (int i=0;i<mx*my;i++){
665 if(bathy->GetValue(i,0)==0){
666 m1->SetValue(i,2,m1col->GetValue(i,0));
667 }
668 else{
669 m1->SetValue(i,2,m0->GetValue(i,2));
670 }
671 if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
672 m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
673 }
674 }
675 delete m1gr;
676 delete m1grsm;
677 delete m1col;
678}/*}}}*/
679double signe(double a){/*{{{*/
680 if(a<0){return -1;}
681 else{return 1;}
682}/*}}}*/
683void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my){/*{{{*/
[18439]684 A->MatrixEqual(Ain);
[18435]685 for (int i=1;i<my-1;i++){
686 for(int j=1;j<mx-1;j++){
687 A->SetValue(i,j,(ctr*Ain->GetValue(i,j)+sd*(Ain->GetValue(i-1,j)+Ain->GetValue(i+1,j)+Ain->GetValue(i,j-1)+Ain->GetValue(i,j+1)))/(ctr+4*sd));
688 }
689 }
690}/*}}}*/
[18446]691double coolshed(double T0,int k,double c,double D){/*{{{*/
692 double T1=T0*exp(-c*pow(k,1/D));
693 return T1;
694}/*}}}*/
Note: See TracBrowser for help on using the repository browser.