快速扫描法（Fast Sweeping Method，FSM）及最短路径辅助的快速扫描法SPAFSM走时计算

已于 2024-06-21 09:14:32 修改
阅读量702 收藏 2
点赞数 3
分类专栏： 地球物理 计算物理 计算数学 文章标签： 算法 c++ c语言
于 2023-10-18 19:57:54 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_42262783/article/details/133913451
版权

概述：

代码包含快速扫描法FSM走时计算的基本方法（Zhao 2005）以及高阶FSM方法和对FSM方法的一点改进：添加最短路径辅助解提高计算精度SPAFSM（Zhang 2023）。该方法及其改进方法，仅供初学者参考，也很容易向高维和各向异性走时计算推广，如有疑问，欢迎随时交流讨论。

(PDF) A shortest‐path‐aided fast‐sweeping method to improve the accuracy of traveltime calculation in vertically transverse isotropic media VTI介质最短路径辅助的快速扫描法提高走时计算精度 (researchgate.net)

(PDF) A shortest-path-aided fast sweeping method to improve the accuracy of travel time calculation (researchgate.net)

(PDF) A shortest-path-aided fast sweeping method to improve the accuracy of travel time calculation (researchgate.net)

1. 源码：

/****************************************************************************************/
/***********                          Fast Sweeping Method                               **********/
/***************************Written By Zhang Jianming,2020.10.10****************/
/***************************************CopyRight************************************/

#include<stdio.h>
#include<iostream>
#include<fstream>
#include<math.h>
#include<stdlib.h>
#include<malloc.h>
#include<string.h>
using namespace std;
float pi=3.1415926;
int nx = 101;
int nz = 101;
float dx = 10;
float dz = 10;
void BubbleSort(float* h, size_t len);
void Swap(float& a, float& b);
float T_control(float a)
{
    if(a>0) return  a;
    else return 0;

}
float min(float a,float b)
{
    if (a>b) return b;
    else return a;
}

void main()
{
    int i,j;
    int x0 = nx/2; int z0 = nz/2;
    //int x0 = 20; int z0 = 20;
    //int x0 = 0; int z0 = 0;
    //int x0 = 375; int z0 = 0;
    float **s=new float*[nx],**T_old=new float*[nx],**T_new=new float*[nx],**T=new float*[nx],**delta=new float*[nx],**eta=new float*[nx],**v=new float*[nx],**vt=new float*[nx];
    for(int i=0;i<nx;i++)
    {
        s[i]=new float[nz];T_old[i]=new float[nz];T_new[i]=new float[nz];T[i]=new float[nz];delta[i]=new float[nz];eta[i]=new float[nz];v[i]=new float[nz];vt[i]=new float[nz];
        for(int j=0;j<nz;j++)
        {    s[i][j]=0;T_old[i][j]=0;T_new[i][j]=0;T[i][j]=0;delta[i][j]=0;eta[i][j]=0;v[i][j]=0;vt[i][j]=0;}
    }
    


    char name[256];
    FILE *fp;
    FILE *filename;
    //filename=fopen("mar_140_501_h15s","rb");
    //for(i=0;i<nx;i++)
    //    for(j=0;j<nz;j++)
    //    fread(&v[i][j], sizeof(float), 1, filename);
    //fclose(filename);
    
    for(i=0;i<nx;i++)   
        for(j=0;j<nz;j++)   
        {
            v[i][j]=1000.0;
            //v[i][j]=1000.0+1*sqrt((i-x0)*dx*(i-x0)*dx+(j-z0)*dz*(j-z0)*dz);
            //if(j>60) v[i][j]=2000.0;
            s[i][j]=(1.0/v[i][j]);
            T_old[i][j]=10.0;
            T[i][j]=22.0;
            delta[i][j]=0.0;
            vt[i][j]=0.0;
            eta[i][j]=0.0;
        }  
    T_old[x0][z0]=0.0;
    //T_old[0][0]=0.0;
    //T_old[x0-1][z0]=dx*s[x0][z0];
    //T_old[x0+1][z0]=dx*s[x0][z0];
    //T_old[x0][z0-1]=dx*s[x0][z0];
    //T_old[x0][z0+1]=dx*s[x0][z0];
    //T_old[x0-1][z0-1]=sqrt(2.0)*dx*s[x0][z0];
    //T_old[x0+1][z0+1]=sqrt(2.0)*dx*s[x0][z0];
    //T_old[x0+1][z0-1]=sqrt(2.0)*dx*s[x0][z0];
    //T_old[x0-1][z0+1]=sqrt(2.0)*dx*s[x0][z0];

    float T_xmin=33.0; float T_zmin=33.0;float A=0.0;float B=0.0;float C=0.0,D=0.0;
    float tmp;
    int loop=0,Max_loop=25;
    int sx,sz;
    int xa,xb,za,zb;
    int xa1,xb1,za1,zb1;

    //SGFSM, notice: it is not an orthogonal coordinate system when dx!=dz
    float tt1=clock();
    for(loop=0;loop<Max_loop;loop++)
    {
        for(i=0;i<nx;i++)
            for(j=0;j<nz;j++)
            {
                xa=i+1;xb=i-1;za=j+1;zb=j-1;
                if(xa==nx) xa=i;
                if(xb==-1) xb=i;
                if(za==nz) za=j;
                if(zb==-1) zb=j;

                1stFD
                //T_xmin=min(T_old[xa][j],T_old[xb][j]);
                //T_zmin=min(T_old[i][za],T_old[i][zb]);
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx;
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                else T[i][j]=(sqrt(2.0)/2.0*T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
                
                //T_xmin=min(T_old[xa][za],T_old[xb][zb]);
                //T_zmin=min(T_old[xb][za],T_old[xa][zb]);
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx*sqrt(2.0);
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*sqrt(2.0)*dx*s[i][j]*sqrt(2.0)*dx-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
            

                //float a[3]={T_xmin,T_zmin,min(T_xmin,T_zmin)};
                //BubbleSort(a, 3);
                //float W=a[0];float V=a[1];float U=a[2];
                //T[i][j]=W+s[i][j]*dx;
                //if(T[i][j]>V) 
                //{
                //    T[i][j]=((V+W+sqrt(2.0*pow(s[i][j]*dx,2)-pow(V-W,2))))/2.0;
                //}
                //if(T[i][j]>U)
                //{
                //    T[i][j]=(V+W+U+sqrt(3.0*pow(s[i][j]*dx,2)-pow(V-W,2)-pow(V-U,2)-pow(U-W,2)))/3.0;
                //}
                //T_old[i][j]=min(T[i][j],T_old[i][j]);

                //T_old[i][j]=min(T_old[xb][zb]+s[i][j]*dx*sqrt(2),T_old[i][j]);
                //xa1=i+2;xb1=i-2;za1=j+2;zb1=j-2;
                //if(xa1==nx) xa1=xa;
                //if(xb1==-1) xb1=xb;
                //if(za1==nz) za1=za;
                //if(zb1==-1) zb1=zb;
                //if(xa1==nx+1) xa1=i;
                //if(xb1==-2) xb1=i;
                //if(za1==nz+1) za1=j;
                //if(zb1==-2) zb1=j;
                //T_old[i][j]=min(T_old[xb1][zb]+s[i][j]*dx*sqrt(5),T_old[i][j]);
                //T_old[i][j]=min(T_old[xb][zb1]+s[i][j]*dx*sqrt(5),T_old[i][j]);


                if(T_old[i][j]>=T_old[xb][zb]&&T_old[i][j]>=T_old[i][zb]&&T_old[i][j]>=T_old[xb][j])
                //if(i>=x0&&j>=z0)
                {
                    //T_old[i][zb]=min(T_old[xb][zb]+s[i][j]*dx,T_old[i][zb]);
                    //T_old[xb][j]=min(T_old[xb][zb]+s[i][j]*dx,T_old[xb][j]);
                    //T_old[i][zb]=min(T_old[xb][zb]+s[i][zb]*dx,T_old[i][zb]);
                    //T_old[xb][j]=min(T_old[xb][zb]+s[xb][j]*dx,T_old[xb][j]);
                    //T_old[i][zb]=min(T_old[xb][zb]+s[xb][zb]*dx,T_old[i][zb]);
                    //T_old[xb][j]=min(T_old[xb][zb]+s[xb][zb]*dx,T_old[xb][j]);
                    T_old[i][zb]=min(T_old[xb][zb]+0.5*(s[xb][zb]+s[i][zb])*dx,T_old[i][zb]);
                    T_old[xb][j]=min(T_old[xb][zb]+0.5*(s[xb][zb]+s[xb][j])*dx,T_old[xb][j]);
                    //tmp=pow(s[i][j]*dx*sqrt(2),2)-pow(T_old[i][zb]-T_old[xb][j],2);
                    //tmp=pow(s[xb][zb]*dx*sqrt(2),2)-pow(T_old[i][zb]-T_old[xb][j],2);
                    tmp=pow(0.5*(s[i][j]+s[xb][zb])*dx*sqrt(2),2)-pow(T_old[i][zb]-T_old[xb][j],2);
                    //tmp=pow(0.25*(s[i][j]+s[xb][zb]+s[i][zb]+s[xb][j])*dx*sqrt(2),2)-pow(T_old[i][zb]-T_old[xb][j],2);
                    if(tmp>=0) T_old[i][j]=min(T_old[xb][zb]+sqrt(tmp),T_old[i][j]);
                    
                    //T_old[i][j]=min(T_old[xb][zb]+s[xb][za]*dx*sqrt(2),T_old[i][j]);
                    //tmp=2*s[i][j]*dx*s[i][j]*dx-pow(T_old[i][zb]-T_old[xb][j],2);
                    //if(tmp>0) T_old[i][j]=min((T_old[xb][j]+T_old[i][zb]+sqrt(tmp))/2,T_old[i][j]);
                }

            }

        for(i=nx-1;i>=0;i--)
            for(j=0;j<nz;j++)
            {
                xa=i+1;xb=i-1;za=j+1;zb=j-1;
                if(xa==nx) xa=i;
                if(xb==-1) xb=i;
                if(za==nz) za=j;
                if(zb==-1) zb=j;
                //T_xmin=min(T_old[xa][j],T_old[xb][j]);
                //T_zmin=min(T_old[i][za],T_old[i][zb]);
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx;
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                else T[i][j]=(sqrt(2.0)/2.0*T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
                
                //T_xmin=min(T_old[xa][za],T_old[xb][zb]);
                //T_zmin=min(T_old[xb][za],T_old[xa][zb]);
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx*sqrt(2.0);
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*sqrt(2.0)*dx*s[i][j]*sqrt(2.0)*dx-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
                
                //T_old[i][j]=min(T_old[xa][zb]+s[i][j]*dx*sqrt(2),T_old[i][j]);
                //xa1=i+2;xb1=i-2;za1=j+2;zb1=j-2;
                //if(xa1==nx) xa1=xa;
                //if(xb1==-1) xb1=xb;
                //if(za1==nz) za1=za;
                //if(zb1==-1) zb1=zb;
                //if(xa1==nx+1) xa1=i;
                //if(xb1==-2) xb1=i;
                //if(za1==nz+1) za1=j;
                //if(zb1==-2) zb1=j;
                //T_old[i][j]=min(T_old[xa1][zb]+s[i][j]*dx*sqrt(5),T_old[i][j]);
                //T_old[i][j]=min(T_old[xa][zb1]+s[i][j]*dx*sqrt(5),T_old[i][j]);

                if(T_old[i][j]>=T_old[xa][zb]&&T_old[i][j]>=T_old[xa][j]&&T_old[i][j]>=T_old[i][zb])
                //if(i<=x0&&j>=z0)
                {
                    //T_old[xa][j]=min(T_old[xa][zb]+s[i][j]*dx,T_old[xa][j]);
                    //T_old[i][zb]=min(T_old[xa][zb]+s[i][j]*dx,T_old[i][zb]);
                    //T_old[xa][j]=min(T_old[xa][zb]+s[xa][j]*dx,T_old[xa][j]);
                    //T_old[i][zb]=min(T_old[xa][zb]+s[i][zb]*dx,T_old[i][zb]);
                    //T_old[xa][j]=min(T_old[xa][zb]+s[xa][zb]*dx,T_old[xa][j]);
                    //T_old[i][zb]=min(T_old[xa][zb]+s[xa][zb]*dx,T_old[i][zb]);
                    T_old[xa][j]=min(T_old[xa][zb]+0.5*(s[xa][zb]+s[xa][j])*dx,T_old[xa][j]);
                    T_old[i][zb]=min(T_old[xa][zb]+0.5*(s[xa][zb]+s[i][zb])*dx,T_old[i][zb]);
                    //tmp=pow(s[i][j]*dx*sqrt(2),2)-pow(T_old[i][zb]-T_old[xa][j],2);
                    //tmp=pow(s[xa][zb]*dx*sqrt(2),2)-pow(T_old[i][zb]-T_old[xa][j],2);
                    tmp=pow(0.5*(s[i][j]+s[xa][zb])*dx*sqrt(2),2)-pow(T_old[i][zb]-T_old[xa][j],2);
                    //tmp=pow(0.25*(s[i][j]+s[xa][zb]+s[xa][j]+s[i][zb])*dx*sqrt(2),2)-pow(T_old[i][zb]-T_old[xa][j],2);
                    if(tmp>=0) T_old[i][j]=min(T_old[xa][zb]+sqrt(tmp),T_old[i][j]);
                    
                    //T_old[i][j]=min(T_old[xa][zb]+s[xb][za]*dx*sqrt(2),T_old[i][j]);
                    //tmp=2*s[i][j]*dx*s[i][j]*dx-pow(T_old[i][zb]-T_old[xa][j],2);
                    //if(tmp>0) T_old[i][j]=min((T_old[xa][j]+T_old[i][zb]+sqrt(tmp))/2,T_old[i][j]);
                }
            }


        for(i=nx-1;i>=0;i--)
            for(j=nz-1;j>=0;j--)
            {
                xa=i+1;xb=i-1;za=j+1;zb=j-1;
                if(xa==nx) xa=i;
                if(xb==-1) xb=i;
                if(za==nz) za=j;
                if(zb==-1) zb=j;
                //T_xmin=min(T_old[xa][j],T_old[xb][j]);
                //T_zmin=min(T_old[i][za],T_old[i][zb]);
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx;
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                else T[i][j]=(sqrt(2.0)/2.0*T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
                
                //T_xmin=min(T_old[xa][za],T_old[xb][zb]);
                //T_zmin=min(T_old[xb][za],T_old[xa][zb]);
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx*sqrt(2.0);
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*sqrt(2.0)*dx*s[i][j]*sqrt(2.0)*dx-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
                
                //T_old[i][j]=min(T_old[xa][za]+s[i][j]*dx*sqrt(2),T_old[i][j]);
                //xa1=i+2;xb1=i-2;za1=j+2;zb1=j-2;
                //if(xa1==nx) xa1=xa;
                //if(xb1==-1) xb1=xb;
                //if(za1==nz) za1=za;
                //if(zb1==-1) zb1=zb;
                //if(xa1==nx+1) xa1=i;
                //if(xb1==-2) xb1=i;
                //if(za1==nz+1) za1=j;
                //if(zb1==-2) zb1=j;
                //T_old[i][j]=min(T_old[xa1][za]+s[i][j]*dx*sqrt(5),T_old[i][j]);
                //T_old[i][j]=min(T_old[xa][za1]+s[i][j]*dx*sqrt(5),T_old[i][j]);

                if(T_old[i][j]>=T_old[xa][za]&&T_old[i][j]>=T_old[xa][j]&&T_old[i][j]>=T_old[i][za])
                //if(i<=x0&&j<=z0)
                {
                    //T_old[xa][j]=min(T_old[xa][za]+s[i][j]*dx,T_old[xa][j]);
                    //T_old[i][za]=min(T_old[xa][za]+s[i][j]*dx,T_old[i][za]);
                    //T_old[xa][j]=min(T_old[xa][za]+s[xa][j]*dx,T_old[xa][j]);
                    //T_old[i][za]=min(T_old[xa][za]+s[i][za]*dx,T_old[i][za]);
                    //T_old[xa][j]=min(T_old[xa][za]+s[xa][za]*dx,T_old[xa][j]);
                    //T_old[i][za]=min(T_old[xa][za]+s[xa][za]*dx,T_old[i][za]);
                    T_old[xa][j]=min(T_old[xa][za]+0.5*(s[xa][za]+s[xa][j])*dx,T_old[xa][j]);
                    T_old[i][za]=min(T_old[xa][za]+0.5*(s[xa][za]+s[i][za])*dx,T_old[i][za]);
                    //tmp=pow(s[i][j]*dx*sqrt(2),2)-pow(T_old[i][za]-T_old[xa][j],2);
                    //tmp=pow(s[xa][za]*dx*sqrt(2),2)-pow(T_old[i][za]-T_old[xa][j],2);
                    tmp=pow(0.5*(s[i][j]+s[xa][za])*dx*sqrt(2),2)-pow(T_old[i][za]-T_old[xa][j],2);
                    //tmp=pow(0.25*(s[i][j]+s[xa][za]+s[xa][j]+s[i][za])*dx*sqrt(2),2)-pow(T_old[i][za]-T_old[xa][j],2);
                    if(tmp>=0) T_old[i][j]=min(T_old[xa][za]+sqrt(tmp),T_old[i][j]);
                    
                    //T_old[i][j]=min(T_old[xa][za]+s[xb][za]*dx*sqrt(2),T_old[i][j]);
                    //tmp=2*s[i][j]*dx*s[i][j]*dx-pow(T_old[i][za]-T_old[xa][j],2);
                    //if(tmp>0) T_old[i][j]=min((T_old[xa][j]+T_old[i][za]+sqrt(tmp))/2,T_old[i][j]);
                }
            }

        for(i=0;i<nx;i++)
            for(j=nz-1;j>=0;j--)
            {
                xa=i+1;xb=i-1;za=j+1;zb=j-1;
                if(xa==nx) xa=i;
                if(xb==-1) xb=i;
                if(za==nz) za=j;
                if(zb==-1) zb=j;
                //T_xmin=min(T_old[xa][j],T_old[xb][j]);
                //T_zmin=min(T_old[i][za],T_old[i][zb]);
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx;
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                else T[i][j]=(sqrt(2.0)/2.0*T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
                
                //T_xmin=min(T_old[xa][za],T_old[xb][zb]);
                //T_zmin=min(T_old[xb][za],T_old[xa][zb]);
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx*sqrt(2.0);
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*sqrt(2.0)*dx*s[i][j]*sqrt(2.0)*dx-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
                
                //T_old[i][j]=min(T_old[xb][za]+s[i][j]*dx*sqrt(2),T_old[i][j]);
                //xa1=i+2;xb1=i-2;za1=j+2;zb1=j-2;
                //if(xa1==nx) xa1=xa;
                //if(xb1==-1) xb1=xb;
                //if(za1==nz) za1=za;
                //if(zb1==-1) zb1=zb;
                //if(xa1==nx+1) xa1=i;
                //if(xb1==-2) xb1=i;
                //if(za1==nz+1) za1=j;
                //if(zb1==-2) zb1=j;
                //T_old[i][j]=min(T_old[xb1][za]+s[i][j]*dx*sqrt(5),T_old[i][j]);
                //T_old[i][j]=min(T_old[xb][za1]+s[i][j]*dx*sqrt(5),T_old[i][j]);
            
                if(T_old[i][j]>=T_old[xb][za]&&T_old[i][j]>=T_old[xb][j]&&T_old[i][j]>=T_old[i][za])
                //if(i>=x0&&j<=z0)
                {
                    //T_old[xb][j]=min(T_old[xb][za]+s[i][j]*dx,T_old[xb][j]);
                    //T_old[i][za]=min(T_old[xb][za]+s[i][j]*dx,T_old[i][za]);
                    //T_old[xb][j]=min(T_old[xb][za]+s[xb][j]*dx,T_old[xb][j]);
                    //T_old[i][za]=min(T_old[xb][za]+s[i][za]*dx,T_old[i][za]);
                    //T_old[xb][j]=min(T_old[xb][za]+s[xb][za]*dx,T_old[xb][j]);
                    //T_old[i][za]=min(T_old[xb][za]+s[xb][za]*dx,T_old[i][za]);
                    T_old[xb][j]=min(T_old[xb][za]+0.5*(s[xb][za]+s[xb][j])*dx,T_old[xb][j]);
                    T_old[i][za]=min(T_old[xb][za]+0.5*(s[xb][za]+s[i][za])*dx,T_old[i][za]);
                    //tmp=pow(s[i][j]*dx*sqrt(2),2)-pow(T_old[i][za]-T_old[xb][j],2);
                    //tmp=pow(s[xb][za]*dx*sqrt(2),2)-pow(T_old[i][za]-T_old[xb][j],2);
                    tmp=pow(0.5*(s[i][j]+s[xb][za])*dx*sqrt(2),2)-pow(T_old[i][za]-T_old[xb][j],2);
                    //tmp=pow(0.25*(s[i][j]+s[xb][za]+s[xb][j]+s[i][za])*dx*sqrt(2),2)-pow(T_old[i][za]-T_old[xb][j],2);
                    if(tmp>=0) T_old[i][j]=min(T_old[xb][za]+sqrt(tmp),T_old[i][j]);
                    
                    //T_old[i][j]=min(T_old[xb][za]+s[xb][za]*dx*sqrt(2),T_old[i][j]);
                    //tmp=2*s[i][j]*dx*s[i][j]*dx-pow(T_old[i][za]-T_old[xb][j],2);
                    //if(tmp>0) T_old[i][j]=min((T_old[xb][j]+T_old[i][za]+sqrt(tmp))/2,T_old[i][j]);
                }
            }
    }//end of loop of 1st FSM 
    float tt2=clock();
    cout<<"One shot modeling for 1stFSM is "<<(double)(tt2-tt1)/CLOCKS_PER_SEC<<" s."<<endl;
    

#if 0
    //2nd order FSM
    float Txf;float Txz;float Tzf; float Tzz;
    float wxf,wxz,wzf,wzz,rxf,rxz,rzf,rzz; float epsi=1e-6;
    //int xa1,xb1,za1,zb1;
    for(loop=0;loop<Max_loop;loop++)
    //for(loop=0;loop<25;loop++)
    {
        for(i=0;i<nx;i++)
            for(j=0;j<nz;j++)
            {
                xa=i+1;xb=i-1;za=j+1;zb=j-1;
                if(xa==nx) xa=i;
                if(xb==-1) xb=i;
                if(za==nz) za=j;
                if(zb==-1) zb=j;
                
                xa1=i+2;xb1=i-2;za1=j+2;zb1=j-2;
                if(xa1==nx) xa1=xa;
                if(xb1==-1) xb1=xb;
                if(za1==nz) za1=za;
                if(zb1==-1) zb1=zb;
                if(xa1==nx+1) xa1=i;
                if(xb1==-2) xb1=i;
                if(za1==nz+1) za1=j;
                if(zb1==-2) zb1=j;
            
                rxf=(epsi+pow(T_old[i][j]-2*T_old[xb][j]+T_old[xb1][j],2))/(epsi+pow(T_old[xa][j]-2*T_old[i][j]+T_old[xb][j],2));
                rxz=(epsi+pow(T_old[xa1][j]-2*T_old[xa][j]+T_old[i][j],2))/(epsi+pow(T_old[xa][j]-2*T_old[i][j]+T_old[xb][j],2));
                rzf=(epsi+pow(T_old[i][j]-2*T_old[i][zb]+T_old[i][zb1],2))/(epsi+pow(T_old[i][za]-2*T_old[i][j]+T_old[i][zb],2));
                rzz=(epsi+pow(T_old[i][za1]-2*T_old[i][za]+T_old[i][j],2))/(epsi+pow(T_old[i][za]-2*T_old[i][j]+T_old[i][zb],2));
                wxf=1.0/(1.0+2.0*pow(rxf,2));
                wxz=1.0/(1.0+2.0*pow(rxz,2));
                wzf=1.0/(1.0+2.0*pow(rzf,2));
                wzz=1.0/(1.0+2.0*pow(rzz,2));
                Txf=(1.0-wxf)*(T_old[xa][j]-T_old[xb][j])/2.0/dx+wxf*(3.0*T_old[i][j]-4.0*T_old[xb][j]+T_old[xb1][j])/2.0/dx;
                Txz=(1.0-wxz)*(T_old[xa][j]-T_old[xb][j])/2.0/dx+wxz*(-3.0*T_old[i][j]+4.0*T_old[xa][j]-T_old[xa1][j])/2.0/dx;
                Tzf=(1.0-wzf)*(T_old[i][za]-T_old[i][zb])/2.0/dz+wzf*(3.0*T_old[i][j]-4.0*T_old[i][zb]+T_old[i][zb1])/2.0/dz;
                Tzz=(1.0-wzz)*(T_old[i][za]-T_old[i][zb])/2.0/dz+wzz*(-3.0*T_old[i][j]+4.0*T_old[i][za]-T_old[i][za1])/2.0/dz;
                T_xmin=min(T_old[i][j]-dx*Txf,T_old[i][j]+dx*Txz);
                T_zmin=min(T_old[i][j]-dz*Tzf,T_old[i][j]+dz*Tzz);
                if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx;
                else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                T_old[i][j]=min(T[i][j],T_old[i][j]);
                
                //rxf=(epsi+pow(T_old[i][j]-2*T_old[xb][zb]+T_old[xb1][zb1],2))/(epsi+pow(T_old[xa][za]-2*T_old[i][j]+T_old[xb][zb],2));
                //wxf=1.0/(1.0+2.0*pow(rxf,2));
                //Txf=(1.0-wxf)*(T_old[xa][za]-T_old[xb][zb])/2.0/sqrt(2)/dx+wxf*(3.0*T_old[i][j]-4.0*T_old[xb][zb]+T_old[xb1][zb1])/2.0/sqrt(2)/dx;
                //T_old[i][j]=min(T_old[i][j]-sqrt(2)*dx*Txf+s[i][j]*dx*sqrt(2),T_old[i][j]);
                
                //rxf=(epsi+pow(T_old[i][j]-2*T_old[xb][zb]+T_old[xb1][zb1],2))/(epsi+pow(T_old[xa][za]-2*T_old[i][j]+T_old[xb][zb],2));
                //rxz=(epsi+pow(T_old[xa1][za1]-2*T_old[xa][za]+T_old[i][j],2))/(epsi+pow(T_old[xa][za]-2*T_old[i][j]+T_old[xb][zb],2));
                //rzf=(epsi+pow(T_old[i][j]-2*T_old[xa][zb]+T_old[xa1][zb1],2))/(epsi+pow(T_old[xb][za]-2*T_old[i][j]+T_old[xa][zb],2));
                //rzz=(epsi+pow(T_old[xb1][za1]-2*T_old[xb][za]+T_old[i][j],2))/(epsi+pow(T_old[xb][za]-2*T_old[i][j]+T_old[xa][zb],2));
                //wxf=1.0/(1.0+2.0*pow(rxf,2));
                //wxz=1.0/(1.0+2.0*pow(rxz,2));
                //wzf=1.0/(1.0+2.0*pow(rzf,2));
                //wzz=1.0/(1.0+2.0*pow(rzz,2));
                //Txf=(1.0-wxf)*(T_old[xa][za]-T_old[xb][zb])/2.0/sqrt(2)/dx+wxf*(3.0*T_old[i][j]-4.0*T_old[xb][zb]+T_old[xb1][zb1])/2.0/sqrt(2)/dx;
                //Txz=(1.0-wxz)*(T_old[xa][za]-T_old[xb][zb])/2.0/sqrt(2)/dx+wxz*(-3.0*T_old[i][j]+4.0*T_old[xa][za]-T_old[xa1][za1])/2.0/sqrt(2)/dx;
                //Tzf=(1.0-wzf)*(T_old[xb][za]-T_old[xa][zb])/2.0/sqrt(2)/dz+wzf*(3.0*T_old[i][j]-4.0*T_old[xa][zb]+T_old[xa1][zb1])/2.0/sqrt(2)/dz;
                //Tzz=(1.0-wzz)*(T_old[xb][za]-T_old[xa][zb])/2.0/sqrt(2)/dz+wzz*(-3.0*T_old[i][j]+4.0*T_old[xb][za]-T_old[xb1][za1])/2.0/sqrt(2)/dz;
                //T_xmin=min(T_old[i][j]-dx*sqrt(2)*Txf,T_old[i][j]+dx*sqrt(2)*Txz);
                //T_zmin=min(T_old[i][j]-dz*sqrt(2)*Tzf,T_old[i][j]+dz*sqrt(2)*Tzz);
                cout<<"i="<<i<<",  j="<<j<<"  "<<Txf<<"  "<<Txz<<"  "<<Tzf<<"  "<<Tzz<<endl;
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx*sqrt(2)) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx*sqrt(2);
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*sqrt(2)*s[i][j]*dx*sqrt(2)-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
            }

        //#if 0
        for(i=nx-1;i>=0;i--)
            for(j=0;j<nz;j++)
            {
                xa=i+1;xb=i-1;za=j+1;zb=j-1;
                if(xa==nx) xa=i;
                if(xb==-1) xb=i;
                if(za==nz) za=j;
                if(zb==-1) zb=j;
                
                xa1=i+2;xb1=i-2;za1=j+2;zb1=j-2;
                if(xa1==nx) xa1=xa;
                if(xb1==-1) xb1=xb;
                if(za1==nz) za1=za;
                if(zb1==-1) zb1=zb;
                if(xa1==nx+1) xa1=i;
                if(xb1==-2) xb1=i;
                if(za1==nz+1) za1=j;
                if(zb1==-2) zb1=j;
            
                rxf=(epsi+pow(T_old[i][j]-2*T_old[xb][j]+T_old[xb1][j],2))/(epsi+pow(T_old[xa][j]-2*T_old[i][j]+T_old[xb][j],2));
                rxz=(epsi+pow(T_old[xa1][j]-2*T_old[xa][j]+T_old[i][j],2))/(epsi+pow(T_old[xa][j]-2*T_old[i][j]+T_old[xb][j],2));
                rzf=(epsi+pow(T_old[i][j]-2*T_old[i][zb]+T_old[i][zb1],2))/(epsi+pow(T_old[i][za]-2*T_old[i][j]+T_old[i][zb],2));
                rzz=(epsi+pow(T_old[i][za1]-2*T_old[i][za]+T_old[i][j],2))/(epsi+pow(T_old[i][za]-2*T_old[i][j]+T_old[i][zb],2));
                wxf=1.0/(1.0+2.0*pow(rxf,2));
                wxz=1.0/(1.0+2.0*pow(rxz,2));
                wzf=1.0/(1.0+2.0*pow(rzf,2));
                wzz=1.0/(1.0+2.0*pow(rzz,2));
                Txf=(1.0-wxf)*(T_old[xa][j]-T_old[xb][j])/2.0/dx+wxf*(3.0*T_old[i][j]-4.0*T_old[xb][j]+T_old[xb1][j])/2.0/dx;
                Txz=(1.0-wxz)*(T_old[xa][j]-T_old[xb][j])/2.0/dx+wxz*(-3.0*T_old[i][j]+4.0*T_old[xa][j]-T_old[xa1][j])/2.0/dx;
                Tzf=(1.0-wzf)*(T_old[i][za]-T_old[i][zb])/2.0/dz+wzf*(3.0*T_old[i][j]-4.0*T_old[i][zb]+T_old[i][zb1])/2.0/dz;
                Tzz=(1.0-wzz)*(T_old[i][za]-T_old[i][zb])/2.0/dz+wzz*(-3.0*T_old[i][j]+4.0*T_old[i][za]-T_old[i][za1])/2.0/dz;
                T_xmin=min(T_old[i][j]-dx*Txf,T_old[i][j]+dx*Txz);
                T_zmin=min(T_old[i][j]-dz*Tzf,T_old[i][j]+dz*Tzz);
                if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx;
                else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                T_old[i][j]=min(T[i][j],T_old[i][j]);
                
                //rxf=(epsi+pow(T_old[i][j]-2*T_old[xb][zb]+T_old[xb1][zb1],2))/(epsi+pow(T_old[xa][za]-2*T_old[i][j]+T_old[xb][zb],2));
                //rxz=(epsi+pow(T_old[xa1][za1]-2*T_old[xa][za]+T_old[i][j],2))/(epsi+pow(T_old[xa][za]-2*T_old[i][j]+T_old[xb][zb],2));
                //rzf=(epsi+pow(T_old[i][j]-2*T_old[xa][zb]+T_old[xa1][zb1],2))/(epsi+pow(T_old[xb][za]-2*T_old[i][j]+T_old[xa][zb],2));
                //rzz=(epsi+pow(T_old[xb1][za1]-2*T_old[xb][za]+T_old[i][j],2))/(epsi+pow(T_old[xb][za]-2*T_old[i][j]+T_old[xa][zb],2));
                //wxf=1.0/(1.0+2.0*pow(rxf,2));
                //wxz=1.0/(1.0+2.0*pow(rxz,2));
                //wzf=1.0/(1.0+2.0*pow(rzf,2));
                //wzz=1.0/(1.0+2.0*pow(rzz,2));
                //Txf=(1.0-wxf)*(T_old[xa][za]-T_old[xb][zb])/2.0/sqrt(2)/dx+wxf*(3.0*T_old[i][j]-4.0*T_old[xb][zb]+T_old[xb1][zb1])/2.0/sqrt(2)/dx;
                //Txz=(1.0-wxz)*(T_old[xa][za]-T_old[xb][zb])/2.0/sqrt(2)/dx+wxz*(-3.0*T_old[i][j]+4.0*T_old[xa][za]-T_old[xa1][za1])/2.0/sqrt(2)/dx;
                //Tzf=(1.0-wzf)*(T_old[xb][za]-T_old[xa][zb])/2.0/sqrt(2)/dz+wzf*(3.0*T_old[i][j]-4.0*T_old[xa][zb]+T_old[xa1][zb1])/2.0/sqrt(2)/dz;
                //Tzz=(1.0-wzz)*(T_old[xb][za]-T_old[xa][zb])/2.0/sqrt(2)/dz+wzz*(-3.0*T_old[i][j]+4.0*T_old[xb][za]-T_old[xb1][za1])/2.0/sqrt(2)/dz;
                //T_xmin=min(T_old[i][j]-dx*sqrt(2)*Txf,T_old[i][j]+dx*sqrt(2)*Txz);
                //T_zmin=min(T_old[i][j]-dz*sqrt(2)*Tzf,T_old[i][j]+dz*sqrt(2)*Tzz);
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx*sqrt(2)) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx*sqrt(2);
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*sqrt(2)*s[i][j]*dx*sqrt(2)-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
            }
        
        for(i=nx-1;i>=0;i--)
            for(j=nz-1;j>=0;j--)
            {
                xa=i+1;xb=i-1;za=j+1;zb=j-1;
                if(xa==nx) xa=i;
                if(xb==-1) xb=i;
                if(za==nz) za=j;
                if(zb==-1) zb=j;
                
                xa1=i+2;xb1=i-2;za1=j+2;zb1=j-2;
                if(xa1==nx) xa1=xa;
                if(xb1==-1) xb1=xb;
                if(za1==nz) za1=za;
                if(zb1==-1) zb1=zb;
                if(xa1==nx+1) xa1=i;
                if(xb1==-2) xb1=i;
                if(za1==nz+1) za1=j;
                if(zb1==-2) zb1=j;
            
                rxf=(epsi+pow(T_old[i][j]-2*T_old[xb][j]+T_old[xb1][j],2))/(epsi+pow(T_old[xa][j]-2*T_old[i][j]+T_old[xb][j],2));
                rxz=(epsi+pow(T_old[xa1][j]-2*T_old[xa][j]+T_old[i][j],2))/(epsi+pow(T_old[xa][j]-2*T_old[i][j]+T_old[xb][j],2));
                rzf=(epsi+pow(T_old[i][j]-2*T_old[i][zb]+T_old[i][zb1],2))/(epsi+pow(T_old[i][za]-2*T_old[i][j]+T_old[i][zb],2));
                rzz=(epsi+pow(T_old[i][za1]-2*T_old[i][za]+T_old[i][j],2))/(epsi+pow(T_old[i][za]-2*T_old[i][j]+T_old[i][zb],2));
                wxf=1.0/(1.0+2.0*pow(rxf,2));
                wxz=1.0/(1.0+2.0*pow(rxz,2));
                wzf=1.0/(1.0+2.0*pow(rzf,2));
                wzz=1.0/(1.0+2.0*pow(rzz,2));
                Txf=(1.0-wxf)*(T_old[xa][j]-T_old[xb][j])/2.0/dx+wxf*(3.0*T_old[i][j]-4.0*T_old[xb][j]+T_old[xb1][j])/2.0/dx;
                Txz=(1.0-wxz)*(T_old[xa][j]-T_old[xb][j])/2.0/dx+wxz*(-3.0*T_old[i][j]+4.0*T_old[xa][j]-T_old[xa1][j])/2.0/dx;
                Tzf=(1.0-wzf)*(T_old[i][za]-T_old[i][zb])/2.0/dz+wzf*(3.0*T_old[i][j]-4.0*T_old[i][zb]+T_old[i][zb1])/2.0/dz;
                Tzz=(1.0-wzz)*(T_old[i][za]-T_old[i][zb])/2.0/dz+wzz*(-3.0*T_old[i][j]+4.0*T_old[i][za]-T_old[i][za1])/2.0/dz;
                T_xmin=min(T_old[i][j]-dx*Txf,T_old[i][j]+dx*Txz);
                T_zmin=min(T_old[i][j]-dz*Tzf,T_old[i][j]+dz*Tzz);
                if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx;
                else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                T_old[i][j]=min(T[i][j],T_old[i][j]);
                
                //rxf=(epsi+pow(T_old[i][j]-2*T_old[xb][zb]+T_old[xb1][zb1],2))/(epsi+pow(T_old[xa][za]-2*T_old[i][j]+T_old[xb][zb],2));
                //rxz=(epsi+pow(T_old[xa1][za1]-2*T_old[xa][za]+T_old[i][j],2))/(epsi+pow(T_old[xa][za]-2*T_old[i][j]+T_old[xb][zb],2));
                //rzf=(epsi+pow(T_old[i][j]-2*T_old[xa][zb]+T_old[xa1][zb1],2))/(epsi+pow(T_old[xb][za]-2*T_old[i][j]+T_old[xa][zb],2));
                //rzz=(epsi+pow(T_old[xb1][za1]-2*T_old[xb][za]+T_old[i][j],2))/(epsi+pow(T_old[xb][za]-2*T_old[i][j]+T_old[xa][zb],2));
                //wxf=1.0/(1.0+2.0*pow(rxf,2));
                //wxz=1.0/(1.0+2.0*pow(rxz,2));
                //wzf=1.0/(1.0+2.0*pow(rzf,2));
                //wzz=1.0/(1.0+2.0*pow(rzz,2));
                //Txf=(1.0-wxf)*(T_old[xa][za]-T_old[xb][zb])/2.0/sqrt(2)/dx+wxf*(3.0*T_old[i][j]-4.0*T_old[xb][zb]+T_old[xb1][zb1])/2.0/sqrt(2)/dx;
                //Txz=(1.0-wxz)*(T_old[xa][za]-T_old[xb][zb])/2.0/sqrt(2)/dx+wxz*(-3.0*T_old[i][j]+4.0*T_old[xa][za]-T_old[xa1][za1])/2.0/sqrt(2)/dx;
                //Tzf=(1.0-wzf)*(T_old[xb][za]-T_old[xa][zb])/2.0/sqrt(2)/dz+wzf*(3.0*T_old[i][j]-4.0*T_old[xa][zb]+T_old[xa1][zb1])/2.0/sqrt(2)/dz;
                //Tzz=(1.0-wzz)*(T_old[xb][za]-T_old[xa][zb])/2.0/sqrt(2)/dz+wzz*(-3.0*T_old[i][j]+4.0*T_old[xb][za]-T_old[xb1][za1])/2.0/sqrt(2)/dz;
                //T_xmin=min(T_old[i][j]-dx*sqrt(2)*Txf,T_old[i][j]+dx*sqrt(2)*Txz);
                //T_zmin=min(T_old[i][j]-dz*sqrt(2)*Tzf,T_old[i][j]+dz*sqrt(2)*Tzz);
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx*sqrt(2)) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx*sqrt(2);
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*sqrt(2)*s[i][j]*dx*sqrt(2)-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
            }

        for(i=0;i<nx;i++)
            for(j=nz-1;j>=0;j--)
            {
                xa=i+1;xb=i-1;za=j+1;zb=j-1;
                if(xa==nx) xa=i;
                if(xb==-1) xb=i;
                if(za==nz) za=j;
                if(zb==-1) zb=j;
                
                xa1=i+2;xb1=i-2;za1=j+2;zb1=j-2;
                if(xa1==nx) xa1=xa;
                if(xb1==-1) xb1=xb;
                if(za1==nz) za1=za;
                if(zb1==-1) zb1=zb;
                if(xa1==nx+1) xa1=i;
                if(xb1==-2) xb1=i;
                if(za1==nz+1) za1=j;
                if(zb1==-2) zb1=j;
            
                rxf=(epsi+pow(T_old[i][j]-2*T_old[xb][j]+T_old[xb1][j],2))/(epsi+pow(T_old[xa][j]-2*T_old[i][j]+T_old[xb][j],2));
                rxz=(epsi+pow(T_old[xa1][j]-2*T_old[xa][j]+T_old[i][j],2))/(epsi+pow(T_old[xa][j]-2*T_old[i][j]+T_old[xb][j],2));
                rzf=(epsi+pow(T_old[i][j]-2*T_old[i][zb]+T_old[i][zb1],2))/(epsi+pow(T_old[i][za]-2*T_old[i][j]+T_old[i][zb],2));
                rzz=(epsi+pow(T_old[i][za1]-2*T_old[i][za]+T_old[i][j],2))/(epsi+pow(T_old[i][za]-2*T_old[i][j]+T_old[i][zb],2));
                wxf=1.0/(1.0+2.0*pow(rxf,2));
                wxz=1.0/(1.0+2.0*pow(rxz,2));
                wzf=1.0/(1.0+2.0*pow(rzf,2));
                wzz=1.0/(1.0+2.0*pow(rzz,2));
                Txf=(1.0-wxf)*(T_old[xa][j]-T_old[xb][j])/2.0/dx+wxf*(3.0*T_old[i][j]-4.0*T_old[xb][j]+T_old[xb1][j])/2.0/dx;
                Txz=(1.0-wxz)*(T_old[xa][j]-T_old[xb][j])/2.0/dx+wxz*(-3.0*T_old[i][j]+4.0*T_old[xa][j]-T_old[xa1][j])/2.0/dx;
                Tzf=(1.0-wzf)*(T_old[i][za]-T_old[i][zb])/2.0/dz+wzf*(3.0*T_old[i][j]-4.0*T_old[i][zb]+T_old[i][zb1])/2.0/dz;
                Tzz=(1.0-wzz)*(T_old[i][za]-T_old[i][zb])/2.0/dz+wzz*(-3.0*T_old[i][j]+4.0*T_old[i][za]-T_old[i][za1])/2.0/dz;
                T_xmin=min(T_old[i][j]-dx*Txf,T_old[i][j]+dx*Txz);
                T_zmin=min(T_old[i][j]-dz*Tzf,T_old[i][j]+dz*Tzz);
                if (fabs(T_xmin-T_zmin)>=s[i][j]*dx) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx;
                else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*s[i][j]*dx-pow(T_xmin-T_zmin,2)))/2;
                T_old[i][j]=min(T[i][j],T_old[i][j]);
                
                //rxf=(epsi+pow(T_old[i][j]-2*T_old[xb][zb]+T_old[xb1][zb1],2))/(epsi+pow(T_old[xa][za]-2*T_old[i][j]+T_old[xb][zb],2));
                //rxz=(epsi+pow(T_old[xa1][za1]-2*T_old[xa][za]+T_old[i][j],2))/(epsi+pow(T_old[xa][za]-2*T_old[i][j]+T_old[xb][zb],2));
                //rzf=(epsi+pow(T_old[i][j]-2*T_old[xa][zb]+T_old[xa1][zb1],2))/(epsi+pow(T_old[xb][za]-2*T_old[i][j]+T_old[xa][zb],2));
                //rzz=(epsi+pow(T_old[xb1][za1]-2*T_old[xb][za]+T_old[i][j],2))/(epsi+pow(T_old[xb][za]-2*T_old[i][j]+T_old[xa][zb],2));
                //wxf=1.0/(1.0+2.0*pow(rxf,2));
                //wxz=1.0/(1.0+2.0*pow(rxz,2));
                //wzf=1.0/(1.0+2.0*pow(rzf,2));
                //wzz=1.0/(1.0+2.0*pow(rzz,2));
                //Txf=(1.0-wxf)*(T_old[xa][za]-T_old[xb][zb])/2.0/sqrt(2)/dx+wxf*(3.0*T_old[i][j]-4.0*T_old[xb][zb]+T_old[xb1][zb1])/2.0/sqrt(2)/dx;
                //Txz=(1.0-wxz)*(T_old[xa][za]-T_old[xb][zb])/2.0/sqrt(2)/dx+wxz*(-3.0*T_old[i][j]+4.0*T_old[xa][za]-T_old[xa1][za1])/2.0/sqrt(2)/dx;
                //Tzf=(1.0-wzf)*(T_old[xb][za]-T_old[xa][zb])/2.0/sqrt(2)/dz+wzf*(3.0*T_old[i][j]-4.0*T_old[xa][zb]+T_old[xa1][zb1])/2.0/sqrt(2)/dz;
                //Tzz=(1.0-wzz)*(T_old[xb][za]-T_old[xa][zb])/2.0/sqrt(2)/dz+wzz*(-3.0*T_old[i][j]+4.0*T_old[xb][za]-T_old[xb1][za1])/2.0/sqrt(2)/dz;
                //T_xmin=min(T_old[i][j]-dx*sqrt(2)*Txf,T_old[i][j]+dx*sqrt(2)*Txz);
                //T_zmin=min(T_old[i][j]-dz*sqrt(2)*Tzf,T_old[i][j]+dz*sqrt(2)*Tzz);
                //if (fabs(T_xmin-T_zmin)>=s[i][j]*dx*sqrt(2)) T[i][j]=min(T_xmin,T_zmin)+s[i][j]*dx*sqrt(2);
                //else T[i][j]=(T_xmin+T_zmin+sqrt(2*s[i][j]*dx*sqrt(2)*s[i][j]*dx*sqrt(2)-pow(T_xmin-T_zmin,2)))/2;
                //T_old[i][j]=min(T[i][j],T_old[i][j]);
            }
    }//end of 2nd order FSM
    float tt3=clock();
    cout<<"One shot modeling for 2ndFSM is "<<(double)(tt3-tt1)/CLOCKS_PER_SEC<<" s."<<endl;
#endif


    //sprintf(name,"V.dat");
    //sprintf(name,"traveltime_cgv_SGFSM.dat");
    sprintf(name,"traveltime_SGFSM.dats");
    //sprintf(name,"traveltime_stagered_FSM.dat");
    //sprintf(name,"traveltime_FSM.dat");
    //sprintf(name,"traveltime_True.dat");
    //fp=fopen("traveltime.dat","wb+");
    //sprintf(name,"traveltime_Mar_stagered_FSM.dat");
    fp=fopen(name,"wb+");
    for(i=0;i<nx;i++)
        for(j=0;j<nz;j++)
        {
            //T_new[i][j]=dx*sqrt(pow((i-x0),2)+pow((j-z0),2))*s[i][j];
            float z=(1.0+1.0/2.0*s[i][j]/1000.0*1.0*((i-x0)*dx*(i-x0)*dx+(j-z0)*dz*(j-z0)*dz));
            T_new[i][j]=log(z+sqrt(z*z-1));
            //T_old[i][j]=T_old[i][j]-T_new[i][j];
            fwrite(&T_old[i][j],sizeof(float),1,fp);
            //fwrite(&T_new[i][j],sizeof(float),1,fp);
            //fwrite(&v[i][j],sizeof(float),1,fp);
        }
    fclose(fp);

    for(i=0;i<nx;i++)
    {
        delete [] s[i];delete [] T_old[i];delete [] T_new[i];delete [] T[i];delete [] delta[i];delete [] eta[i];delete [] v[i];delete [] vt[i];
    }
    delete [] s;delete [] T_old;delete [] T_new;delete [] T;delete [] delta;delete [] eta;delete [] v;delete [] vt;

}

void BubbleSort(float* h, size_t len)
{
    if(h==NULL) return;
    if(len<=1) return;
    //i是次数，j是具体下标 
    for(int i=0;i<len-1;++i)
        for(int j=0;j<len-1-i;++j)
            if(h[j]>h[j+1])
                Swap(h[j],h[j+1]);               
                return;
} 
void Swap(float& a, float& b)
{
    float t=a;
    a=b;
    b=t;
    return;
}

2. 数值算例

2.1 均匀速度模型计算的走时场对比：

1阶：

高阶：

2.2 复杂模型（Marmousi）的算例：

张建明0704
关注
专栏目录
三维重建之平面扫描算法（Plane-sweeping）
程序员，他们想的是什么?他们想的永远都是技术，他们崇尚的也永远都是技术。
01-25 248
三维重建之平面扫描算法（Plane-sweeping）
C#快速扫描线算法——求取倾斜四边形与图片像素点交叠区域
Eco-dream‘s CODE World
10-30 2373
1.问题简要说明： a.图片的像素点大小为1×1的小矩形，以图片左下角为坐标原点，建立图片的直角坐标系 b.四边形的左下点在Y轴上半轴任意位置，宽为Width(>100)，高为Height(>10),且与X轴有负小角度夹角( 问题： 求取四边形所影响的像素点，及影响大小（四边形完全覆盖像素点时，定义影响大小为1） 2.快速扫描算法 a.传统扫描算法（C#） 1）根据四边形左下顶点的坐
fast sweeping method
wavefront24的博客
11-22 1534
对于一类H-J方程，可运用FSM进行求解 扫描顺序：  for(sx=-1;sx      for(sy=-1;sy         for(i=(sx=0:i             for(j=(sy=0:j 对所有网格节点上的函数值φ，进行a,b,c三个步骤的搜索，在搜索过程中，当si>0时，每个节点函数值的搜索顺序是从小到大；si   算法具体思路： %所有节点的φ(i,
各向异性TTI介质快速扫描法（Fast Sweeping Method，FSM）走时计算
最新发布
weixin_42262783的博客
11-09 421
本人做了复现，虽然是复现，但是也花了很久的时间。其中的一些问题，想了很久才解决。当然实现了以后，就觉得所有的问题原来都不是什么特别的难题。2.沿x轴和z轴群速度的计算，可以采用Shooting方法（Qian，2012），或者Fomel（2004）近似公式。1.因果律的物理含义，代表群速度（射线方向）要和走时更新的方向一致，其实代表了射线的传播方向，即群速度方向。上面这俩式子很重要，又能用来推导沿x轴和z轴群速度，还能用来做因果律的判定。直接上结果了，代码有需要的私聊我。
快速扫描算法提取鱼眼图像有效区域
wd1603926823的专栏
05-12 4870
这个原理还是来自于《鱼眼图像校正算法研究》 这是作者在这篇论文中提出的自己改进的算法（我姑且叫它快速扫描算法吧） 我觉得这个算法中有个地方写错了 可能是作者笔误吧 第二步：p（x，y）与p（x+1，y）的亮度之差如果大于T，则p是左边界的切点 我觉得这里不是这样 我按照这个算法编程 实现出来 效果不对 我觉得第二步应该是p（x+1，y）大于T，则p（x，y）是左边界的切点。我实现出来效果才正确。因
计算几何中的扫描（sweeping)算法（转自百度）
CandyLien
09-05 2704
要解决的问题是判断一大堆线段中是否存在一对线段相交。 个人感觉计算几何的扫描算法算是属于解题思想上的东西，类似常规问题的搜索，是遍历处理对象的一种方法。 算法的思想是安照某种方式组织要被处理的对象，比如点、线段。组织方式通常就是从上到下、从左到右等。具体的办法是假想一
快速行进法
chh0018的博客
12-29 3288
快速行进法 Eikonal 方程（程函方程）的数值求解方法有 快速行进法(Fast Marching Method) 快速扫描法(Fast Sweeping Method)，适用范围更广，效率更高。 快速行进法(Fast Marching Method） 简介 快速行进算法 （Fast Marching Method） 是求解程函方程 （Eikonal Equation）F ∣ ∇ T ∣ = 1 的一种高效数值算法，程函方程属于非线性偏微分方程，可以认为是一种近似波动方程。 在界面演化问题中的求解的
Fast processing method to generate gigabyte computer generated holography for three-dimensional dynamic holographic display
02-13
The computation of CGH for a three-dimensional (3D) dynamic holographic display can also be sped up by programming with our proposed method. It can optimize both file loading and the inline ...
电子功用-基于络合作用与在线扫集技术的毛细管电泳直接紫外法快速检测Amadori化合物
09-15
《电子功用-基于络合作用与在线扫集技术的毛细管电泳直接紫外法快速检测Amadori化合物》 本文主要探讨了在现代生物分析领域中，如何利用电子技术来提升毛细管电泳（Capillary Electrophoresis, CE）的性能，特别是...
横扫式聚类和扫描法聚类各有什么特点
04-26
横扫式聚类（Sweeping Clustering）是一种快速的聚类算法，它按照数据点在某一方向（如水平方向）上的坐标值进行排序，然后将相邻的数据点归为一组，形成聚类。该算法速度快，但对于不同方向的聚类效果并不理想。 ...
使用Julia进行图像处理--图像分割
盼小辉丶的博客
12-10 4649
使用Julia进行图像处理--图像分割。使用Julia的ImageSegmentation包，学习如何使用有监督和无监督的方法来简化或表示图像，使其更易于分析。有时，图像分割被称为图像识别中的第一步。
Eikonal方程
weixin_33737134的博客
10-31 2218
今天写的内容是关于在level set方法中用到的Eikonal方程。具体讲解如何求解这个方程。 方程的定义：其中d(x)是指x到S的距离 简单介绍一下fast sweeping method。它采用了差分方法和轮换次序的Gauss-Seidel迭代来解决这个离散问题的。我们今天的描述在2D中进行。这种算法分为3步。 我们将需要求解的区域进行网格化，使用xi,j代表计算区域中的网格点，h...
操作系统 磁盘调度---扫描（SCAN）算法 （电梯调度算法）（C++实现 操作系统实验）
weixin_57202646的博客
05-25 3982
思路： 不仅考虑到欲访问的磁道与当前磁道间的距离，更优先考虑的是磁头当前的移动方向，当磁头正在自里向外移动时，SCAN算法所考虑的下一个访问对象应是该磁头自里向外方向上距离最近的那个，直至扫描完最外面的磁头，然后方向变为自外向里扫描 代码如下： #include<iostream> #include<cstring> #include<algorithm> #include<queue> #include<cmath> using n
IEEE Transactions on Antennas and Propagation 202107论文速递
sukeee的博客
07-03 4636
1. An Atomic Receiver for AM and FM Radio Communication 用于 AM 和 FM 无线电通信的原子接收器 David Alexander Anderson. Rachel Elizabeth Sapiro. Georg Raithel. 原文链接 “ Radio reception relies on antennas for the collection of electromagnetic fields carrying informati..
基于python的三维射线追踪库-ttcrpy详解（4）
Neverlevsun的博客
05-02 2421
基于python的三维射线追踪库-ttcrpy详解（4） python库ttcrpy中有三种射线追踪算法，他们是：FSM（快速扫描算法）、SPM（最短路径算法）、DSPM（动节点最短路径算法）。本内容为对比这三种算法的特点。 文章目录基于python的三维射线追踪库-ttcrpy详解（4）1、模型一：低速异常体1、FSM（快速扫描算法）2、SPM（最短路径法）3、DSPM（动节点最短路径算法）4、三种方法追踪得到的走时对比5、模型二、高速异常体6、SPM与DSPM计算时间对比7、python代码搬砖不易，走
求最短路径的四种方法（Dijkstra，Floyd，Bellman-Ford，SPFA算法）
菜鸡成长史
01-28 5278
【前言】 还不知道图是什么的可以看看这篇：图的基本概念。 最短路径，即在网络（带权的图）中，求两个不同顶点之间的所有路径中，边的权值之和最小的那一条路径。第一个顶点为源点，最后一个顶点为终点。 【问题分类】 单源最短路问题（特定点到所有点的最短路）和多源最短路问题（任意两点之间的最短路）。其中又分为有向图和无向图，有权图和无权图。根据边权的正负，又分为带负权边和不带负权边的最短路。 【最...
快速排序算法实现：单向扫描法&&双向扫描法&&三分法(有相同元素)
qq_42837890的博客
03-06 964
单向指针扫描 void Swap(int i, int j, int array[])//交换数组内两个元素 { int temp; temp = array[i]; array[i] = array[j]; array[j] = temp; } void QuickSort(int begin, int end, int array[])//快速排序 { if (beg...
最短寻道时间优先算法（SSTF）&&扫描算法（SCAN）
热门推荐
理科男同学
10-19 5万+
最短寻道时间优先算法（SSTF） SSTF问题描述：SSTF算法选择调度处理的磁道是与当前磁头所在磁道距离最近的磁道，以使每次的寻找时间最短。当然，总是选择最小寻找时间并不能保证平均寻找时间最小，但是能提供比FCFS算法更好的性能。这种算法会产生“饥饿”现象。 1、算法思想：优先选择距当前磁头最近的访问请求进行服务，主要考虑寻道优先。 2、优点：改善了磁盘平均服务时间。 3、缺点：...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值