readhess.c 2.88 KB
Newer Older
Simone Butzer's avatar
Simone Butzer committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
/*------------------------------------------------------------------------
 * Copyright (C) 2015 For the list of authors, see file AUTHORS.
 *
 * This file is part of IFOS3D.
 * 
 * IFOS3D is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
 * 
 * IFOS3D is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with IFOS3D. See file COPYING and/or 
 * <http://www.gnu.org/licenses/gpl-2.0.html>.
--------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
 *   Read Hessian from files  
 *  ----------------------------------------------------------------------*/

#include "fd.h"

26
void readhess(int nx, int ny, int nz, float ***  hess1, float ***  hess2, float ***hess3, float finv, int iteration){
Simone Butzer's avatar
Simone Butzer committed
27 28

	extern int NX, NY, NZ, NXG, NYG, NZG, POS[4];
29
	extern float VP0, VS0, RHO0;
Simone Butzer's avatar
Simone Butzer committed
30
	extern FILE *FP;
31
	extern char  HESS_FILE[STRING_SIZE];
Simone Butzer's avatar
Simone Butzer committed
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
		
	/* local variables */
	float rhov, vp, vs; /*qp, qs;*/
	int i, j, k, ii, jj, kk;
	FILE *fp_vs, *fp_vp, *fp_rho;
	char filename[STRING_SIZE];

	/* choose data format: ASCII: format=2, BINARY: format=3*/
	const int format=3; 
	

		

	   fprintf(FP,"\n...reading hess information from hess-files...\n");

47
	   /*sprintf(filename,"hess/hess.vp");*/
48
	   sprintf(filename,"%s.vp_%4.2fHz_it%d",HESS_FILE,finv,iteration);
Simone Butzer's avatar
Simone Butzer committed
49
	   fp_vp=fopen(filename,"r");
50
	   if (fp_vp==NULL) err(" Could not open hess_vp ! ");
Simone Butzer's avatar
Simone Butzer committed
51

52
	   sprintf(filename,"%s.vs_%4.2fHz_it%d",HESS_FILE,finv,iteration);
Simone Butzer's avatar
Simone Butzer committed
53
	   fp_vs=fopen(filename,"r");
54
	   if (fp_vs==NULL) err(" Could not open  hess_vs! ");
Simone Butzer's avatar
Simone Butzer committed
55

56
	   sprintf(filename,"%s.rho_%4.2fHz_it%d",HESS_FILE,finv,iteration);
Simone Butzer's avatar
Simone Butzer committed
57
	   fp_rho=fopen(filename,"r");
58
	   if (fp_rho==NULL) err(" Could not open hess_rho ! ");
Simone Butzer's avatar
Simone Butzer committed
59 60
	   

61
	
Simone Butzer's avatar
Simone Butzer committed
62
	/* loop over global grid */
63
		for (k=1;k<=NZG;k++){
Simone Butzer's avatar
Simone Butzer committed
64
			for (i=1;i<=NXG;i++){
65 66
				for (j=1;j<=NYG;j++){

Simone Butzer's avatar
Simone Butzer committed
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
				vp=readdsk(fp_vp, format);
				vs=readdsk(fp_vs, format);
				rhov=readdsk(fp_rho , format);
	
				/* only the PE which belongs to the current global gridpoint 
				is saving model parameters in his local arrays */
				if ((POS[1]==((i-1)/NX)) && 
				    (POS[2]==((j-1)/NY)) && 
				    (POS[3]==((k-1)/NZ))){
					ii=i-POS[1]*NX;
					jj=j-POS[2]*NY;
					kk=k-POS[3]*NZ;
					
					/*hess1[jj][ii][kk]=vp;
					hess2[jj][ii][kk]=vs;
					hess3[jj][ii][kk]=rhov;*/
					
					
85 86 87
					hess1[jj][ii][kk]=vp*pow(VP0,2);
					hess2[jj][ii][kk]=vs*pow(VS0,2);
					hess3[jj][ii][kk]=rhov*pow(RHO0,2);
Simone Butzer's avatar
Simone Butzer committed
88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
					
			
				}
				}
			}
		}

	fclose(fp_vp);
	fclose(fp_vs);
	fclose(fp_rho);	
	

}