Changeset 18435

Timestamp:

08/18/14 16:21:52 (11 years ago)

Author:

abuzzi

Message:

MODIF: generation of the new matrix

File:

• issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp (modified) (3 diffs)

issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp

@@ -131 +131 @@
 double 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);
 void GSLsquarefit(Matrix** pX,Matrix* A,Matrix* B);
+double signe(double a);
+void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my);
 /*}}}*/
 
@@ -281 +283 @@
         M->Echo();
         o->Echo();*/
-
+double a;
+a=signe(-9);
+cout<<a;
         return 0;
 }/*}}}*/
@@ -491 +495 @@
         return e;
 }/*}}}*/
+void 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){/*{{{*/
+        Matrix* m1gr=new Matrix(my,mx);
+        Matrix* m1grsm=new Matrix(my,mx);
+        Matrix* m1col=new Matrix(mx*my,1);
+        double u=0;
+        double y=0;
+        m1->MatrixEqual(m0);
+        for (int i=0;i<mx*my;i++){
+                if(icethick->GetValue(i,0)==0){
+                        u=double (rand())/ double(RAND_MAX);
+                        y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
+                        m1->SetValue(i,1,m0->GetValue(i,1)+y*ptval);
+                        if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
+                                m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
+                        }
+                        if(m1->GetValue(i,1)>=m1->GetValue(i,0)+mmax){
+                                m1->SetValue(i,1,m1->GetValue(i,0)+mmax);
+                        }
+                }
+        }
+        m1->ExtractColumn(m1col,1);
+        vec2gridsimple(m1col,m1gr,mx,my);
+        filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
+        reshape(m1grsm,m1col,mx,my);
+        for (int i=0;i<mx*my;i++){
+                if(icethick->GetValue(i,0)==0){
+                        m1->SetValue(i,1,m1col->GetValue(i,0));
+                }
+                else{
+                        m1->SetValue(i,1,m0->GetValue(i,1));
+                }
+                if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
+                        m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
+                }
+        }
+
+        for (int i=0;i<mx*my;i++){
+                if(bathy->GetValue(i,0)==0){
+                        u=double (rand())/ double(RAND_MAX);
+                        y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
+                        m1->SetValue(i,2,m0->GetValue(i,2)+y*ptval);
+                        if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
+                                m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
+                        }
+                        if(m1->GetValue(i,2)>=m1->GetValue(i,1)+mmax2){
+                                m1->SetValue(i,2,m1->GetValue(i,1)+mmax2);
+                        }
+                }
+        }
+        m1->ExtractColumn(m1col,2);
+        vec2gridsimple(m1col,m1gr,mx,my);
+        filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
+        reshape(m1grsm,m1col,mx,my);
+        for (int i=0;i<mx*my;i++){
+                if(bathy->GetValue(i,0)==0){
+                        m1->SetValue(i,2,m1col->GetValue(i,0));
+                }
+                else{
+                        m1->SetValue(i,2,m0->GetValue(i,2));
+                }
+                if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
+                        m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
+                }
+        }
+        delete m1gr;
+        delete m1grsm;
+        delete m1col;
+}/*}}}*/
+double signe(double a){/*{{{*/
+        if(a<0){return -1;}
+        else{return 1;}
+}/*}}}*/
+void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my){/*{{{*/
+        for (int i=1;i<my-1;i++){
+                for(int j=1;j<mx-1;j++){
+                        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));
+                }
+        }
+}/*}}}*/