devbox
很楠不爱3
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

三维各向同性快速扫描法3DFSM

作者：很楠不爱3 | 2024-03-26 14:18:43

赞

踩

三维各向同性快速扫描法3DFSM

之前快速扫描法的帖子中提供的二维代码稍显凌乱，为此追加提供一个新的三维的代码，更加清楚直观地展示三维FSM的过程，便于初学者学习

代码：

  1. /****************************************************************************************/
  2. /***********                         3D Fast Sweeping Method                          **********/
  3. /**************************Written By Zhang Jianming, 2020.10.10****************/
  4. /***************************************CopyRight************************************/
  5.  
  6. #include<stdio.h>
  7. #include<iostream>
  8. #include<fstream>
  9. #include<math.h>
  10. #include<stdlib.h>
  11. #include<malloc.h>
  12. #include<string.h>
  13. #include<ctime>
  14. using namespace std;
  15. float pi=3.1415926;
  16. int nx =101;
  17. int ny=101;
  18. int nz =101;
  19. float dx= 10;
  20. float dy= 10;
  21. float dz =10;
  22. int T_max=1e3;
  23. int T_min=0;
  24. float T_control(float a);
  25. float min(float a,float b);
  26. float max(float a,float b);
  27. void BubbleSort(float* h, size_t len);
  28. void Swap(float& a, float& b);
  29.  
  30.  
  31. int main()
  32. {
  33. 	int i,j,k;
  34. 	int x0 = nx/2;int y0=ny/2; int z0 = nz/2;
  35. 	float ***s3=new float** [nx];
  36. 	float ***T_old3=new float** [nx];
  37. 	float ***T3=new float** [nx];
  38. 	float ***v3=new float** [nx];
  39. 	for(int i=0;i<nx;i++)
  40. 	{
  41. 		s3[i]=new float *[ny];
  42. 		T_old3[i]=new float *[ny];
  43. 		T3[i]=new float *[ny];
  44. 		v3[i]=new float *[ny];
  45. 		for(int j=0;j<ny;j++)
  46. 		{
  47. 			s3[i][j]=new float[nz];
  48. 			T_old3[i][j]=new float[nz];
  49. 			v3[i][j]=new float[nz];
  50. 			T3[i][j]=new float[nz];
  51. 		}
  52. 	}
  53.  
  54.  
  55.  
  56. 	//initial for 3D
  57. 	for(i=0;i<nx;i++)   
  58. 		for(j=0;j<ny;j++)
  59. 			for(k=0;k<nz;k++)
  60. 			{
  61. 				v3[i][j][k]=1000;
  62. 				s3[i][j][k]=(1.0/v3[i][j][k]);
  63. 				T3[i][j][k]=22.0;
  64. 				T_old3[i][j][k]=22.0;
  65. 			}
  66. 	int sx=x0;int sy=y0;int sz=z0;
  67. 	T_old3[sx][sy][sz]=0.0;
  68.  
  69. 	float T_xmin=33.0;float T_ymin=33.0; float T_zmin=33.0;float A=0.0;float B=0.0;float C=0.0,D=0.0;
  70. 	int loop=0,Max_loop=1;
  71. 	/***********FP FSM for 3D*****************/
  72. 	float tt1=clock();
  73. 	for(loop=0;loop<Max_loop;loop++)
  74. 	{
  75. 		for(i=sx;i<nx;i++)
  76. 			for(j=sy;j<ny;j++)
  77. 				for(k=sz;k<nz;k++)
  78. 				{
  79. 					int xa=i+1;if(xa==nx) xa=nx-1;
  80. 					int xb=i-1;if(xb==-1) xb=0;
  81. 					int ya=j+1;if(ya==ny) ya=ny-1;
  82. 					int yb=j-1;if(yb==-1) yb=0;
  83. 					int za=k+1;if(za==nz) za=nz-1;
  84. 					int zb=k-1;if(zb==-1) zb=0;
  85. 					T_xmin=min(T_old3[xa][j][k],T_old3[xb][j][k]);
  86. 					T_ymin=min(T_old3[i][ya][k],T_old3[i][yb][k]);
  87. 					T_zmin=min(T_old3[i][j][za],T_old3[i][j][zb]);
  88. 					float a[3]={T_xmin,T_ymin,T_zmin};
  89. 					BubbleSort(a, 3);
  90. 					float W=a[0];float V=a[1];float U=a[2];
  91. 					T3[i][j][k]=W+s3[i][j][k]*dx;
  92. 					if(T3[i][j][k]>V) 
  93. 					{
  94. 						T3[i][j][k]=((V+W+sqrt(2.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2))))/2.0;
  95. 						//T3[i][j][k]=((V+W+sqrt(-pow(V,2)-pow(W,2)+2.0*V*W))/(2.0*s3[i][j][k]*dx*s3[i][j][k]*dx))/2.0;
  96. 					}
  97. 					if(T3[i][j][k]>U)
  98. 					{
  99. 						T3[i][j][k]=(V+W+U+sqrt(3.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2)-pow(V-U,2)-pow(U-W,2)))/3.0;
  100. 						float tssp=T3[xb][yb][zb]+s3[xb][yb][zb]*dx*sqrt(3);
  101. 						T3[i][j][k]=min(T3[i][j][k],tssp);
  102. 						//float b=-2.0*(W+V+U);
  103. 						//T3[i][j][k]=(2.0*(W+V+U)+sqrt(b*b-12*(W*W+V*V+U*U-pow(s3[i][j][k]*dx,2))))/6.0;
  104. 						//T3[i][j][k]=((2.0*pow(V+W+U,2)+sqrt(4.0*pow(V+W+U,2)))/(12.0*pow(W*W+V*V+U*U-s3[i][j][k]*dx,2)))/2.0;
  105. 					}
  106. 					T_old3[i][j][k]=min(T_old3[i][j][k],T3[i][j][k]);
  107. 				}
  108. 		/**/
  109. 		for(i=sx;i>=0;i--)
  110. 			for(j=sy;j<ny;j++)
  111. 				for(k=sz;k<nz;k++)
  112. 				{
  113. 					int xa=i+1;if(xa==nx) xa=nx-1;
  114. 					int xb=i-1;if(xb==-1) xb=0;
  115. 					int ya=j+1;if(ya==ny) ya=ny-1;
  116. 					int yb=j-1;if(yb==-1) yb=0;
  117. 					int za=k+1;if(za==nz) za=nz-1;
  118. 					int zb=k-1;if(zb==-1) zb=0;
  119. 					T_xmin=min(T_old3[xa][j][k],T_old3[xb][j][k]);
  120. 					T_ymin=min(T_old3[i][ya][k],T_old3[i][yb][k]);
  121. 					T_zmin=min(T_old3[i][j][za],T_old3[i][j][zb]);
  122. 					float a[3]={T_xmin,T_ymin,T_zmin};
  123. 					BubbleSort(a, 3);
  124. 					float W=a[0];float V=a[1];float U=a[2];
  125. 					T3[i][j][k]=W+s3[i][j][k]*dx;
  126. 					if(T3[i][j][k]>V) 
  127. 					{
  128. 						T3[i][j][k]=((V+W+sqrt(2.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2))))/2.0;
  129. 						//T3[i][j][k]=((V+W+sqrt(-pow(V,2)-pow(W,2)+2.0*V*W))/(2.0*s3[i][j][k]*dx*s3[i][j][k]*dx))/2.0;
  130. 					}
  131. 					if(T3[i][j][k]>U)
  132. 					{
  133. 						T3[i][j][k]=(V+W+U+sqrt(3.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2)-pow(V-U,2)-pow(U-W,2)))/3.0;
  134. 						float tssp=T3[xa][yb][zb]+s3[xa][yb][zb]*dx*sqrt(3);
  135. 						T3[i][j][k]=min(T3[i][j][k],tssp);
  136. 						//float b=-2.0*(W+V+U);
  137. 						//T3[i][j][k]=(2.0*(W+V+U)+sqrt(b*b-12*(W*W+V*V+U*U-pow(s3[i][j][k]*dx,2))))/6.0;
  138. 						//T3[i][j][k]=((2.0*pow(V+W+U,2)+sqrt(4.0*pow(V+W+U,2)))/(12.0*pow(W*W+V*V+U*U-s3[i][j][k]*dx,2)))/2.0;
  139. 					}
  140. 					T_old3[i][j][k]=min(T_old3[i][j][k],T3[i][j][k]);
  141. 				}
  142.  
  143. 		for(i=sx;i<nx;i++)
  144. 			for(j=sy;j>=0;j--)
  145. 				for(k=sz;k<nz;k++)
  146. 				{
  147. 					int xa=i+1;if(xa==nx) xa=nx-1;
  148. 					int xb=i-1;if(xb==-1) xb=0;
  149. 					int ya=j+1;if(ya==ny) ya=ny-1;
  150. 					int yb=j-1;if(yb==-1) yb=0;
  151. 					int za=k+1;if(za==nz) za=nz-1;
  152. 					int zb=k-1;if(zb==-1) zb=0;
  153. 					T_xmin=min(T_old3[xa][j][k],T_old3[xb][j][k]);
  154. 					T_ymin=min(T_old3[i][ya][k],T_old3[i][yb][k]);
  155. 					T_zmin=min(T_old3[i][j][za],T_old3[i][j][zb]);
  156. 					float a[3]={T_xmin,T_ymin,T_zmin};
  157. 					BubbleSort(a, 3);
  158. 					float W=a[0];float V=a[1];float U=a[2];
  159. 					T3[i][j][k]=W+s3[i][j][k]*dx;
  160. 					if(T3[i][j][k]>V) 
  161. 					{
  162. 						T3[i][j][k]=((V+W+sqrt(2.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2))))/2.0;
  163. 						//T3[i][j][k]=((V+W+sqrt(-pow(V,2)-pow(W,2)+2.0*V*W))/(2.0*s3[i][j][k]*dx*s3[i][j][k]*dx))/2.0;
  164. 					}
  165. 					if(T3[i][j][k]>U)
  166. 					{
  167. 						T3[i][j][k]=(V+W+U+sqrt(3.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2)-pow(V-U,2)-pow(U-W,2)))/3.0;
  168. 						float tssp=T3[xb][ya][zb]+s3[xb][ya][zb]*dx*sqrt(3);
  169. 						T3[i][j][k]=min(T3[i][j][k],tssp);
  170. 						//float b=-2.0*(W+V+U);
  171. 						//T3[i][j][k]=(2.0*(W+V+U)+sqrt(b*b-12*(W*W+V*V+U*U-pow(s3[i][j][k]*dx,2))))/6.0;
  172. 						//T3[i][j][k]=((2.0*pow(V+W+U,2)+sqrt(4.0*pow(V+W+U,2)))/(12.0*pow(W*W+V*V+U*U-s3[i][j][k]*dx,2)))/2.0;
  173. 					}
  174. 					T_old3[i][j][k]=min(T_old3[i][j][k],T3[i][j][k]);
  175. 				}
  176.  
  177. 		for(i=sx;i>=0;i--)
  178. 			for(j=sy;j>=0;j--)
  179. 				for(k=sz;k<nz;k++)
  180. 				{
  181. 					int xa=i+1;if(xa==nx) xa=nx-1;
  182. 					int xb=i-1;if(xb==-1) xb=0;
  183. 					int ya=j+1;if(ya==ny) ya=ny-1;
  184. 					int yb=j-1;if(yb==-1) yb=0;
  185. 					int za=k+1;if(za==nz) za=nz-1;
  186. 					int zb=k-1;if(zb==-1) zb=0;
  187. 					T_xmin=min(T_old3[xa][j][k],T_old3[xb][j][k]);
  188. 					T_ymin=min(T_old3[i][ya][k],T_old3[i][yb][k]);
  189. 					T_zmin=min(T_old3[i][j][za],T_old3[i][j][zb]);
  190. 					float a[3]={T_xmin,T_ymin,T_zmin};
  191. 					BubbleSort(a, 3);
  192. 					float W=a[0];float V=a[1];float U=a[2];
  193. 					T3[i][j][k]=W+s3[i][j][k]*dx;
  194. 					if(T3[i][j][k]>V) 
  195. 					{
  196. 						T3[i][j][k]=((V+W+sqrt(2.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2))))/2.0;
  197. 						//T3[i][j][k]=((V+W+sqrt(-pow(V,2)-pow(W,2)+2.0*V*W))/(2.0*s3[i][j][k]*dx*s3[i][j][k]*dx))/2.0;
  198. 					}
  199. 					if(T3[i][j][k]>U)
  200. 					{
  201. 						T3[i][j][k]=(V+W+U+sqrt(3.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2)-pow(V-U,2)-pow(U-W,2)))/3.0;
  202. 						float tssp=T3[xa][ya][zb]+s3[xa][ya][zb]*dx*sqrt(3);
  203. 						T3[i][j][k]=min(T3[i][j][k],tssp);
  204. 						//float b=-2.0*(W+V+U);
  205. 						//T3[i][j][k]=(2.0*(W+V+U)+sqrt(b*b-12*(W*W+V*V+U*U-pow(s3[i][j][k]*dx,2))))/6.0;
  206. 						//T3[i][j][k]=((2.0*pow(V+W+U,2)+sqrt(4.0*pow(V+W+U,2)))/(12.0*pow(W*W+V*V+U*U-s3[i][j][k]*dx,2)))/2.0;
  207. 					}
  208. 					T_old3[i][j][k]=min(T_old3[i][j][k],T3[i][j][k]);
  209. 				}
  210.  
  211. 		for(i=sx;i<nx;i++)
  212. 			for(j=sy;j<ny;j++)
  213. 				for(k=sz;k>=0;k--)
  214. 				{
  215. 					int xa=i+1;if(xa==nx) xa=nx-1;
  216. 					int xb=i-1;if(xb==-1) xb=0;
  217. 					int ya=j+1;if(ya==ny) ya=ny-1;
  218. 					int yb=j-1;if(yb==-1) yb=0;
  219. 					int za=k+1;if(za==nz) za=nz-1;
  220. 					int zb=k-1;if(zb==-1) zb=0;
  221. 					T_xmin=min(T_old3[xa][j][k],T_old3[xb][j][k]);
  222. 					T_ymin=min(T_old3[i][ya][k],T_old3[i][yb][k]);
  223. 					T_zmin=min(T_old3[i][j][za],T_old3[i][j][zb]);
  224. 					float a[3]={T_xmin,T_ymin,T_zmin};
  225. 					BubbleSort(a, 3);
  226. 					float W=a[0];float V=a[1];float U=a[2];
  227. 					T3[i][j][k]=W+s3[i][j][k]*dx;
  228. 					if(T3[i][j][k]>V) 
  229. 					{
  230. 						T3[i][j][k]=((V+W+sqrt(2.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2))))/2.0;
  231. 						//T3[i][j][k]=((V+W+sqrt(-pow(V,2)-pow(W,2)+2.0*V*W))/(2.0*s3[i][j][k]*dx*s3[i][j][k]*dx))/2.0;
  232. 					}
  233. 					if(T3[i][j][k]>U)
  234. 					{
  235. 						T3[i][j][k]=(V+W+U+sqrt(3.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2)-pow(V-U,2)-pow(U-W,2)))/3.0;
  236. 						float tssp=T3[xb][yb][za]+s3[xb][yb][za]*dx*sqrt(3);
  237. 						T3[i][j][k]=min(T3[i][j][k],tssp);
  238. 						//float b=-2.0*(W+V+U);
  239. 						//T3[i][j][k]=(2.0*(W+V+U)+sqrt(b*b-12*(W*W+V*V+U*U-pow(s3[i][j][k]*dx,2))))/6.0;
  240. 						//T3[i][j][k]=((2.0*pow(V+W+U,2)+sqrt(4.0*pow(V+W+U,2)))/(12.0*pow(W*W+V*V+U*U-s3[i][j][k]*dx,2)))/2.0;
  241. 					}
  242. 					T_old3[i][j][k]=min(T_old3[i][j][k],T3[i][j][k]);
  243. 				}
  244.  
  245. 		for(i=sx;i>=0;i--)
  246. 			for(j=sy;j<ny;j++)
  247. 				for(k=sz;k>=0;k--)
  248. 				{
  249. 					int xa=i+1;if(xa==nx) xa=nx-1;
  250. 					int xb=i-1;if(xb==-1) xb=0;
  251. 					int ya=j+1;if(ya==ny) ya=ny-1;
  252. 					int yb=j-1;if(yb==-1) yb=0;
  253. 					int za=k+1;if(za==nz) za=nz-1;
  254. 					int zb=k-1;if(zb==-1) zb=0;
  255. 					T_xmin=min(T_old3[xa][j][k],T_old3[xb][j][k]);
  256. 					T_ymin=min(T_old3[i][ya][k],T_old3[i][yb][k]);
  257. 					T_zmin=min(T_old3[i][j][za],T_old3[i][j][zb]);
  258. 					float a[3]={T_xmin,T_ymin,T_zmin};
  259. 					BubbleSort(a, 3);
  260. 					float W=a[0];float V=a[1];float U=a[2];
  261. 					T3[i][j][k]=W+s3[i][j][k]*dx;
  262. 					if(T3[i][j][k]>V) 
  263. 					{
  264. 						T3[i][j][k]=((V+W+sqrt(2.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2))))/2.0;
  265. 						//T3[i][j][k]=((V+W+sqrt(-pow(V,2)-pow(W,2)+2.0*V*W))/(2.0*s3[i][j][k]*dx*s3[i][j][k]*dx))/2.0;
  266. 					}
  267. 					if(T3[i][j][k]>U)
  268. 					{
  269. 						T3[i][j][k]=(V+W+U+sqrt(3.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2)-pow(V-U,2)-pow(U-W,2)))/3.0;
  270. 						float tssp=T3[xa][yb][za]+s3[xa][yb][za]*dx*sqrt(3);
  271. 						T3[i][j][k]=min(T3[i][j][k],tssp);
  272. 						//float b=-2.0*(W+V+U);
  273. 						//T3[i][j][k]=(2.0*(W+V+U)+sqrt(b*b-12*(W*W+V*V+U*U-pow(s3[i][j][k]*dx,2))))/6.0;
  274. 						//T3[i][j][k]=((2.0*pow(V+W+U,2)+sqrt(4.0*pow(V+W+U,2)))/(12.0*pow(W*W+V*V+U*U-s3[i][j][k]*dx,2)))/2.0;
  275. 					}
  276. 					T_old3[i][j][k]=min(T_old3[i][j][k],T3[i][j][k]);
  277. 				}
  278.  
  279. 		for(i=sx;i<nx;i++)
  280. 			for(j=sy;j>=0;j--)
  281. 				for(k=sz;k>=0;k--)
  282. 				{
  283. 					int xa=i+1;if(xa==nx) xa=nx-1;
  284. 					int xb=i-1;if(xb==-1) xb=0;
  285. 					int ya=j+1;if(ya==ny) ya=ny-1;
  286. 					int yb=j-1;if(yb==-1) yb=0;
  287. 					int za=k+1;if(za==nz) za=nz-1;
  288. 					int zb=k-1;if(zb==-1) zb=0;
  289. 					T_xmin=min(T_old3[xa][j][k],T_old3[xb][j][k]);
  290. 					T_ymin=min(T_old3[i][ya][k],T_old3[i][yb][k]);
  291. 					T_zmin=min(T_old3[i][j][za],T_old3[i][j][zb]);
  292. 					float a[3]={T_xmin,T_ymin,T_zmin};
  293. 					BubbleSort(a, 3);
  294. 					float W=a[0];float V=a[1];float U=a[2];
  295. 					T3[i][j][k]=W+s3[i][j][k]*dx;
  296. 					if(T3[i][j][k]>V) 
  297. 					{
  298. 						T3[i][j][k]=((V+W+sqrt(2.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2))))/2.0;
  299. 						//T3[i][j][k]=((V+W+sqrt(-pow(V,2)-pow(W,2)+2.0*V*W))/(2.0*s3[i][j][k]*dx*s3[i][j][k]*dx))/2.0;
  300. 					}
  301. 					if(T3[i][j][k]>U)
  302. 					{
  303. 						T3[i][j][k]=(V+W+U+sqrt(3.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2)-pow(V-U,2)-pow(U-W,2)))/3.0;
  304. 						float tssp=T3[xb][ya][za]+s3[xb][ya][za]*dx*sqrt(3);
  305. 						T3[i][j][k]=min(T3[i][j][k],tssp);
  306. 						//float b=-2.0*(W+V+U);
  307. 						//T3[i][j][k]=(2.0*(W+V+U)+sqrt(b*b-12*(W*W+V*V+U*U-pow(s3[i][j][k]*dx,2))))/6.0;
  308. 						//T3[i][j][k]=((2.0*pow(V+W+U,2)+sqrt(4.0*pow(V+W+U,2)))/(12.0*pow(W*W+V*V+U*U-s3[i][j][k]*dx,2)))/2.0;
  309. 					}
  310. 					T_old3[i][j][k]=min(T_old3[i][j][k],T3[i][j][k]);
  311. 				}
  312.  
  313. 		for(i=sx;i>=0;i--)
  314. 			for(j=sy;j>=0;j--)
  315. 				for(k=sz;k>=0;k--)
  316. 				{
  317. 					int xa=i+1;if(xa==nx) xa=nx-1;
  318. 					int xb=i-1;if(xb==-1) xb=0;
  319. 					int ya=j+1;if(ya==ny) ya=ny-1;
  320. 					int yb=j-1;if(yb==-1) yb=0;
  321. 					int za=k+1;if(za==nz) za=nz-1;
  322. 					int zb=k-1;if(zb==-1) zb=0;
  323. 					T_xmin=min(T_old3[xa][j][k],T_old3[xb][j][k]);
  324. 					T_ymin=min(T_old3[i][ya][k],T_old3[i][yb][k]);
  325. 					T_zmin=min(T_old3[i][j][za],T_old3[i][j][zb]);
  326. 					float a[3]={T_xmin,T_ymin,T_zmin};
  327. 					BubbleSort(a, 3);
  328. 					float W=a[0];float V=a[1];float U=a[2];
  329. 					T3[i][j][k]=W+s3[i][j][k]*dx;
  330. 					if(T3[i][j][k]>V) 
  331. 					{
  332. 						T3[i][j][k]=((V+W+sqrt(2.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2))))/2.0;
  333. 						//T3[i][j][k]=((V+W+sqrt(-pow(V,2)-pow(W,2)+2.0*V*W))/(2.0*s3[i][j][k]*dx*s3[i][j][k]*dx))/2.0;
  334. 					}
  335. 					if(T3[i][j][k]>U)
  336. 					{
  337. 						T3[i][j][k]=(V+W+U+sqrt(3.0*pow(s3[i][j][k]*dx,2)-pow(V-W,2)-pow(V-U,2)-pow(U-W,2)))/3.0;
  338. 						float tssp=T3[xa][ya][za]+s3[xa][ya][za]*dx*sqrt(3);
  339. 						T3[i][j][k]=min(T3[i][j][k],tssp);
  340. 						//float b=-2.0*(W+V+U);
  341. 						//T3[i][j][k]=(2.0*(W+V+U)+sqrt(b*b-12*(W*W+V*V+U*U-pow(s3[i][j][k]*dx,2))))/6.0;
  342. 						//T3[i][j][k]=((2.0*pow(V+W+U,2)+sqrt(4.0*pow(V+W+U,2)))/(12.0*pow(W*W+V*V+U*U-s3[i][j][k]*dx,2)))/2.0;
  343. 					}
  344. 					T_old3[i][j][k]=min(T_old3[i][j][k],T3[i][j][k]);
  345. 				}
  346. 	}
  347. 	float tt2=clock();
  348. 	cout<<"One shot modeling for 3DFSM is "<<(double)(tt2-tt1)/CLOCKS_PER_SEC<<" s."<<endl;
  349. 	float *Rdep=new float[nx];
  350.  
  351.  
  352. 	ofstream T3out("T3.dat");
  353. 	for(int i=0;i<nx;i++)
  354. 	{
  355. 		for(int j=0;j<ny;j++)
  356. 		{
  357. 			for(int k=0;k<nz;k++)
  358. 			{
  359. 				//T_old3[i][j][k]=sqrt(pow(i-nx/2,2)+pow(j-ny/2,2)+pow(k-0,2))*dx*1.0/2000.0;
  360. 				T_old3[i][j][k]-=sqrt(pow(i-x0,2)+pow(j-y0,2)+pow(k-z0,2))*dx*1.0/1000.0;
  361. 				T3out.write((char*)&T_old3[i][j][k],sizeof(float));
  362. 			}
  363. 		}
  364. 	}
  365. 	cout<<T_old3[nx-1][ny-1][nz-1]<<endl;;
  366. 	T3out.close();
  367.  
  368. 	for(int i=0;i<nx;i++)
  369. 	{
  370. 		for(int j=0;j<ny;j++)
  371. 		{
  372. 			delete [] s3[i][j];
  373. 			delete [] T_old3[i][j];
  374. 			delete [] v3[i][j];
  375. 			delete [] T3[i][j];
  376. 		}
  377. 		delete []s3[i];
  378. 		delete [] T_old3[i];
  379. 		delete []v3[i];
  380. 		delete [] T3[i];
  381. 	}
  382. 	delete []s3;
  383. 	delete [] T_old3;
  384. 	delete []v3;
  385. 	delete [] T3;
  386. 	return 0;
  387. }//end of main
  388. float T_control(float a)
  389. {
  390. 	if(a>0) return  a;
  391. 	else return 0;
  392.  
  393. }
  394. float min(float a,float b)
  395. {
  396. 	if (a>b) return b;
  397. 	else return a;
  398. }
  399. float max(float a,float b)
  400. {
  401. 	if (a>b) return a;
  402. 	else return b;
  403. }
  404. void BubbleSort(float* h, size_t len)
  405. {
  406. 	if(h==NULL) return;
  407. 	if(len<=1) return;
  408. 	//i是次数，j是具体下标 
  409. 	for(int i=0;i<len-1;++i)
  410. 		for(int j=0;j<len-1-i;++j)
  411. 			if(h[j]>h[j+1])
  412. 				Swap(h[j],h[j+1]);               
  413. 	return;
  414. } 
  415. void Swap(float& a, float& b)
  416. {
  417. 	float t=a;
  418. 	a=b;
  419. 	b=t;
  420. 	return;
  421. }
  422.  
 

计算实例
 

走时场（震源位于上表面的中心位置）
 

 

整体走时误差
 

 

 

下底面走时误差
 


                

        
声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/很楠不爱3/article/detail/317752
推荐阅读
相关标签
Copyright ©  2003-2013 www.wpsshop.cn  版权所有，并保留所有权利。

  
闽ICP备14008679号