adjustable source disk radius; verison 2006b

configure.ac

@@ -1,4 +1,4 @@
-AC_INIT([fomosto-qseis], [2006a])
+AC_INIT([fomosto-qseis], [2006b])
 AM_INIT_AUTOMAKE([-Wall -Werror foreign])
 AC_PROG_F77
 AC_CONFIG_FILES([

doc/examples/qs6inp.dat

@@ -63,11 +63,17 @@
 #    the Nyquist, and so on: the larger this parameter, the smaller the space-
 #    domain aliasing effect, but also the more computation effort); 
 # 4. the factor for suppressing time domain aliasing (> 0 and <= 1) (Note 2).
+# 5. new in version 2006b: radius of the (horizontal) source disk [m].
+#    Using a smaller value leads to less artifacts. This can be very slow if
+#    source and receivers are at same or similar depth. A larger value leads to
+#    quicker convergence. Set to -1.0 to use the automatic setting of earlier
+#    versions of QSEIS.
 #------------------------------------------------------------------------------
  2                                   |int: sw_algorithm;
  0.000  0.000  0.000  0.000          |dble: slw(1-4);
  1.00                                |dble: sample_rate;
  0.01                                |dble: supp_factor;
+ -1.0                                |dble: source_radius;
 #------------------------------------------------------------------------------
 #
 #	        OPTIONS FOR PARTIAL SOLUTIONS

src/Makefile.am

@@ -1,3 +1,3 @@
-bin_PROGRAMS = fomosto_qseis2006a
-fomosto_qseis2006a_SOURCES = qsmatrix0.f bessj0.f bessj1.f bessj.f caxcb.f cdgemp.f cmemcpy.f four1.f qsam2ve.f qsbsj.f qsfftinv.f qsgetinp.f qshksh.f qskern.f qslayer.f qsmain.f qsmultis.f qspsv.f qsqmodel.f qssh.f qssource.f qssublay.f qsve2am.f qswavelet.f qswaveno.f qswvint.f taper.f wavelet.f skip_comments.f qsglobal.h 
+bin_PROGRAMS = fomosto_qseis2006b
+fomosto_qseis2006b_SOURCES = qsmatrix0.f bessj0.f bessj1.f bessj.f caxcb.f cdgemp.f cmemcpy.f four1.f qsam2ve.f qsbsj.f qsfftinv.f qsgetinp.f qshksh.f qskern.f qslayer.f qsmain.f qsmultis.f qspsv.f qsqmodel.f qssh.f qssource.f qssublay.f qsve2am.f qswavelet.f qswaveno.f qswvint.f taper.f wavelet.f skip_comments.f qsglobal.h
 

src/qsgetinp.f

@@ -115,6 +115,9 @@ c
         print *,'factor is replaced by the default value of 1/e.'
       endif
       fi=dlog(suppress)/(2.d0*pi*twindow)
+
+      call skip_comments(unit)
+      read(unit,*)rdisk
 c
 c     partial solution parameters
 c     ===========================

src/qsglobal.h

@@ -142,8 +142,8 @@ c
       common /pathfilter/ pathdepth,iflat,ipath,lpath,isurf,ndtrans
 c
       integer nt,nf
-      double precision dt,df,fi
-      common /sampling/ dt,df,fi,nt,nf
+      double precision dt,df,fi,rdisk
+      common /sampling/ dt,df,fi,nt,nf,rdisk
 c
       double precision tstart,twindow,tau,v0
       common /tparas/ tstart,twindow,tau,v0

src/qswvint.f

@@ -80,7 +80,14 @@ c
       c2dk=dcmplx(2.d0*dk,0.d0)
       cdk2=dcmplx(dk*dk,0.d0)
 c
-      r0=dmax1(cmax/(0.1d0*fcut),1.d0/(200.d0*dk))
+c     radius of horizontal disk
+      if(rdisk.lt.0.d0)then
+c       automatic setting of earlier versions of QSEIS
+        r0=dmax1(cmax/(0.1d0*fcut),1.d0/(200.d0*dk))
+      else
+c       explicit setting from configuration file
+        r0=rdisk
+      endif
 c
       if(fullwave)then
         f=0.d0
@@ -368,4 +375,4 @@ c
       endif
 c
       return
-      end
+      end