Resolve Merge issue

doc/latex/6_Parameter_Definition.tex

@@ -587,7 +587,7 @@ If NORMALIZE is set to 1, the synthetic data and the measured data will be norma
 
 To reduce the memory requirements during an inversion one can define that only every DTINV time sample is used for the calculation of the gradients. To set this parameter appropriately one has to keep in mind the Nyquist criterion to avoid aliasing effects.
 
-
+\newpage
 \section{Step length estimation}
 {\color{blue}{\begin{verbatim}
 "Step length estimation" : "comment", 
@@ -608,8 +608,9 @@ For a better definition of the parabola the improved line search is now trying t
 
 EPS\_SCALE += EPS\_SCALE/SCALEFAC.\\
 
-If a corresponding value epst(3) can be found after STEPMAX forward modellings, DENISE can fit a parabola through the 3 points (L2t(i),epst(i)) and estimates an optimum step length at the minimum of the parabola. If the L2-value L2t(3) after STEPMAX forward models is still smaller than L2t(2) the optimum steplength estimated by parabolic fitting will be not larger than epst(3).
+If a corresponding value epst(3) can be found after STEPMAX forward modellings, DENISE can fit a parabola through the 3 points (L2t(i),epst(i)) and estimates an optimum step length at the minimum of the parabola. If the L2-value L2t(3) after STEPMAX forward models is still smaller than L2t(2) the optimum steplength estimated by parabolic fitting will be not larger than epst(3).\\
 
+Please note: This step length search is only available wenn PCG method (GRAD\_METHOD==1) is used and will be automatically selected. 
 \newpage
 
 

src/Makefile

@@ -165,6 +165,7 @@ DENISE= \
     matcopy_viscac.c \
     prepare_update_p.c \
     readmod_viscac.c \
+    time_window_glob.c \
         filter_frequencies.c
 	
 SNAPMERGE_OBJ = $(SNAPMERGE_SCR:%.c=%.o)

src/apply_workflow.c

@@ -38,6 +38,7 @@ void apply_workflow(float ** workflow,int workflow_lines,int workflow_line_curre
     extern float JOINT_INVERSION_PSV_SH_ALPHA_RHO;
     extern int EPRECOND;
     extern float EPSILON_WE;
+    extern int GRAD_METHOD;
     
     /******************/
     /* Apply Workflow */
@@ -119,7 +120,7 @@ void apply_workflow(float ** workflow,int workflow_lines,int workflow_line_curre
     EPRECOND=workflow[workflow_line_current][11];
     EPSILON_WE=workflow[workflow_line_current][12];
     
-    if(*LBFGS_iter_start==*iter){
+    if(*LBFGS_iter_start==*iter && GRAD_METHOD==2){
         if(MYID==0)printf("\n L-BFGS will be used from iteration %d on.",*LBFGS_iter_start+1);
     }
 }
\ No newline at end of file

src/fd.h

@@ -186,7 +186,7 @@ void  inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int
 
 void  inseis_source_wavelet(float *section, int ns, int ishot);
 
-void  taper(float **sectionpdiff, int ntr, int ns);
+void  taper(float *section, int ns, float fc);
 
 void  output_source_signal(FILE *fp, float **signals, int ns, int seis_form);
 

src/globvar.h

@@ -68,7 +68,8 @@ float VPUPPERLIM, VPLOWERLIM, VSUPPERLIM, VSLOWERLIM, RHOUPPERLIM, RHOLOWERLIM;
 
 float npower, k_max_PML;
 
-int INV_STF, N_STF, N_STF_START;
+int TAPER_STF;
+int INV_STF, N_STF, N_STF_START,TW_IND;
 char PARA[STRING_SIZE];
 
 int TIME_FILT, ORDER, ZERO_PHASE;

src/taper.c

@@ -17,20 +17,21 @@
 -----------------------------------------------------------------------------------------*/
 
 /*------------------------------------------------------------------------
- *   Write seismograms to disk                                  
- *   last update 19/01/02, T. Bohlen
+ *   Taper function now adapted for use in stf.c, only for single traces!
+ *   last update 21/10/15, L. Gassner
  *  ----------------------------------------------------------------------*/
 #include "fd.h"
 #include "segy.h"
 
-void  taper(float **sectionpdiff, int ntr, int ns){
+void  taper(float *section, int ns, float fc){
 
 	/* declaration of extern variables */
-	extern int MYID, TAPERLENGTH;
+	extern int MYID;
+	extern float DT;
 	extern FILE *FP;
 	
 	/* declaration of local variables */
-	int i,j, h;
+	int i,j, h, taperlength, taperduration;
 	int tracl1;
 	float a;
 	float damping, amp;
@@ -39,31 +40,47 @@ void  taper(float **sectionpdiff, int ntr, int ns){
 	window = vector(1,ns);
         amp1 = vector(1,ns);
 	
+	taperlength=(int)(ceil(2.0/fc/DT));
+	taperduration=2*taperlength;
+	
 	/* "Cerjan"-Window */
-        damping=100.0;
+        damping=99.9;
         amp=1.0-damping/100.0;
-        a=sqrt(-log(amp)/((TAPERLENGTH-1)*(TAPERLENGTH-1)));
-        for (i=1;i<=ns;i++){
+	        a=sqrt(-log(amp)/((taperlength-1)*(taperlength-1)));
+        
+	for (i=1;i<=ns;i++){
 		window[i]=1.0;
 		amp1[i]=0.0;
 	}
+	
 	if (MYID==0){
-		fprintf(FP,"\n ntr: %d\n",ntr);
-		fprintf(FP,"\n ns: %d\n",ns);
-		fprintf(FP,"\n TAPERLENGTH: %d\n",TAPERLENGTH);
+		fprintf(FP,"\n fc: %f\n",fc);
+		fprintf(FP,"\n taperlength: %d\n",taperlength);
+	}
+	
+	for (i=1;i<=taperlength;i++){
+		amp1[i]=exp(-(a*a*(taperlength-i)*(taperlength-i)));
+	}
+	
+// 	/* Taper at the beginning of the window*/
+// 	for (i=1;i<=taperlength;i++){
+// 		window[i]=amp1[i];
+// 	}
+	
+	h=1;
+	for (i=taperduration;i>=(taperduration-taperlength+3);i--){
+		window[i]=amp1[h];
+		h++;
 	}
-        for (i=1;i<=TAPERLENGTH;i++){amp1[i]=exp(-(a*a*(TAPERLENGTH-i)*(TAPERLENGTH-i)));}
-        for (i=1;i<=1;i++){window[i]=amp1[i];}
-        h=1;
-        for (i=ns;i>=(ns-TAPERLENGTH+3);i--){window[i]=amp1[h];h++;}
 	
-	/* Apply taper window */
-	for(tracl1=1;tracl1<=ntr;tracl1++){ 
-		
-		for(j=1;j<=ns;j++){
-			sectionpdiff[tracl1][j]*=window[j];
-		}
+	for (i=taperduration;i<=ns;i++){
+		window[i]=amp1[i];
 	}
+	
+	for(j=1;j<=ns;j++){
+		section[j]*=window[j];
+	}
+	
 	free_vector(window,1,ns);
 	free_vector(amp1,1,ns);
 }