!< !Tomo ! !Usage: !Tomo.x out.H ! !INPUT PARAMETERS: ! stdin - sepfile Input velocity model in depth ! ! General params ! verb - logical [.false.] Whether or not to be verbose ! tau - logical [.false.] Whether or not to tau tomography ! only_rays-logical [.false.] Whether or not to only do ray tracing ! !Inversion parameters ! dot_test - logical [.false.] Wheter to run the dot test ! niter_prec - integer [10] Number of itterations of preconditioned tomo ! niter_reg - integer [0] Number of itterations of regularized tomo ! eps - real [1.] Scaling factor ! max_mem - real [60] Maximum memory devoted to holding ray information ! delta - real [.05] How important smoothness outside ray area is ! eps_ray - real [1.] Epsilon on how much to damp model outside rays ! time_errors- sepfile [none] If exists write a file with used time errors ! residual_times- sepfile [none] If exists write a file with resid time errors ! delta_slow - sepfile [none] If exists write out changes in model ! semblance - sepfile File containing depth positioning semblance ! ! ! Reflector ! reflectors - sepfile File containing reflector locations ! (required if pvalues or rays not supplied) ! analyze_ref - integer* [all] number of reflectors to analayze ! !Ray Parameters ! zero_rays - sep3dfile [optional] File containing zero rays (t,cmp,nref) ! px_rays - sep3dfile [optional] File containing offset rays(t,cmp,nref) ! oray_crp - real First cmp to analyze (if rays not provided) ! dray_crp - real Sampling of cmp to analyze (if rays not provided) ! nray_crp - integer Number of cmps to analyze (if rays not provided) ! dt - real [.016] Sampling of raytracing to save ! nt - integer [6./dt+1] Number of raytracing steps to do ! jt - integer [4] How much to subsample saved raytracing results ! npx - integer Number of ray parameters ! opx - real First ray parameters ! dpx - real ray parameter sampling ! ! Tau Parameters ! dtau - real Sampling in tau for the velocity model ! ntau - integer Number of tau samples ! ! !Filter Parameters ! isotropic - logical [.false.] Whether or not use an issotropic smoother ! pfile - sepfile File containg slopes for steering filters ! widthfile - sepfile File containg widths for steering filters ! hwidth - real [3.] Half width of triangle filte ! ampfile - sepfile File containing amplitude for steering filters ! amp - real [.9] Amplitude of steering filter ! smooth - real [2.5] Half width of triangle smoother over ! reflector positions ! smooth_model -real [.true.] Wheter or not to apply smoothing to slowness ! or change in slowness ! ! !Time Errors ! delta_t - sep3dfile File containing travel time errors(hkeys:offset_x,cmp_x,iref) ! !OUTPUT PARAMETERS: ! stdout - sepfile Updated velocity model ! ! ! !Description: ! ! ! !> !------------------------------------------------- ! !Author: Robert Clapp, ESMB 463, 7230253 ! !Date Created:Sun Jul 18 20:12:46 PDT 1999 ! !Purpose: ! ! program Tomo_it use get_tomo_params_mod use handle_rays_mod use reg_mod use tomo_chain_mod use solver_mod use cgstep_mod use tomo_base_mod use dottest use time_error_mod implicit none logical :: verb type(tau_conv) :: tau_pars type(vel) :: vel_pars type(tomo) :: tomo_pars type(reflectors) :: reflector_pars type(sep3d) :: px_rays,zero_rays,errors type(reg) :: reg_pars type(nfilter) :: filt type(semb) :: semb_pars real,pointer :: time_errors(:) real,allocatable :: data(:),model(:),err(:),junk(:),mdiag(:) integer :: nv(2),ierr !number of model points in velocity field real :: dot1(2),dot2(2) integer :: i1,nte,nl,ithread,nthread real,allocatable :: res(:),resm(:) real,pointer :: edge(:) logical :: dot_check !################ !READ PARAMETERS# !################ call initpar() call doc('/homes/sep/bob/papers/thesis/tau/Src/Tomo_2d.rs90') call grab_general_params(verb,tau_pars,tomo_pars) call grab_vel_param(vel_pars) call grab_refl_param(reflector_pars) call grab_errors(semb_pars) call grab_ray_params(tomo_pars,reflector_pars,px_rays,zero_rays) if (.not. tomo_pars%only_rays) then call grab_reg_params(vel_pars,tau_pars,reg_pars) end if !############# !DO PREP WORK# !############# call init_rays(verb) if (.not. tomo_pars%only_rays) then call init_reg(verb,tau_pars) end if call calc_rays_new(vel_pars,reflector_pars,semb_pars,zero_rays,px_rays& &) !calculate rays if (ithread.eq.0) then if (tomo_pars%only_rays) then call exit(0) end if !calculate regularization params call create_reg_filter(filt,reg_pars,tomo_pars,reflector_pars,edge& &) call tomo_base_init(verb,tau_pars,vel_pars,tomo_pars,zero_rays,px_rays& &,nv) !Match rays and errors call create_errors(zero_rays,px_rays,semb_pars,reflector_pars& &,time_errors,vel_pars) call tomo_calc_size() !Calculate our operator blocks allocate(data(size(time_errors)+product(nv)),model(product(nv))& &,stat=ierr) data(1:size(time_errors))=time_errors if (ierr.ne.0) then call erexit("trouble allocating data and model") end if call init_chains(nv,data,tomo_pars%eps,edge) call init_tomo(tomo_pars,data,model,reg_pars) !initialize tomography !where(filt%mis) data(size(time_errors)+1:)=0. !data(size(time_errors)+1:)=0. !initialize scaling call init_hit_count(data,tomo_pars%delta,size(time_errors),size(model& &)) !############# !PERFORM TOMO# !############# if (tomo_pars%dot_test) then write(0,*) "BEGINING DOT TESTS" write(0,*) "tomo_oper" call dot_test(tomo_oper,size(model),size(time_errors),dot1,dot2) write(0,*) "wraper" if (tau_pars%tau) then call dot_test(wraper,size(model),size(model),dot1,dot2) end if write(0,*) "i_op" call dot_test(i_op,size(model),size(model),dot1,dot2) if (tomo_pars%niter_prec>0) then write(0,*) "npoydiv" call dot_test(npolydiv_lop,size(model),size(model),dot1,dot2) write(0,*) "prec_chain" call dot_test(prec_chain,size(model),size(data),dot1,dot2) end if write(0,*) "hit_count" call dot_test(hit_count_oper,size(model),size(model),dot1,dot2) write(0,*) "nhelicon" call dot_test(nhelicon_lop,size(model),size(model),dot1,dot2) if (tau_pars%tau) then write(0,*) "calc_delta_sigma" call dot_test(calc_delta_sigma,size(model),size(model),dot1,dot2& &) end if write(0,*) "scale it" call dot_test(scale_it,size(model),size(model),dot1,dot2) write(0,*) "prec base" call dot_test(prec_base,size(model),size(data),dot1,dot2) end if model=0. if (tomo_pars%resolve_it) then allocate(mdiag(size(model))) call solver(tomo_oper,cgstep,x=model,dat=data(1:size(time_errors& &)),mdiag=mdiag,nort=tomo_pars%nort,niter=tomo_pars%niter_reg& &,verb=.true.) call srite("out",mdiag,size(mdiag)*4) else if (tomo_pars%time_errors) then call write_time_errors("time_errors",time_errors) end if if (tomo_pars%niter_prec>0) then allocate(err(tomo_pars%niter_prec)) allocate(res(size(data)),resm(size(model))) call solver_reg(prec_chain,solv=cgstep,niter=tomo_pars%niter_prec& &,x=model,dat=data,x0=model,err=err,verb=verb,reg=prec_reg& &,nreg=(size(model)),eps=tomo_pars%eps_ray,res=res,resm=resm) ! ierr=hit_count_oper(.false.,.false.,model,resm) call delta_calc(.true.,model,tomo_pars%delta_slow) write(0,*) "ERROR:",err deallocate(err) end if if (tomo_pars%niter_reg>0) then if (tomo_pars%niter_prec>0) then deallocate(res) end if allocate(err(tomo_pars%niter_reg),res(size(data))) call solver_reg(tomo_oper,solv=cgstep,niter=tomo_pars%niter_reg& &,x=model,dat=data(1:size(time_errors)),x0=model,err=err& &,verb=verb,reg=nhelicon_lop,nreg=(size(model)),eps=20.,res=res) ierr=tomo_oper(.false.,.false.,model,time_errors) deallocate(err) end if if (tomo_pars%time_errors) then call write_time_errors("residual_times",res(1:size(time_errors& &))) end if nte=size(time_errors) nl=size(res) i1=nl-nte call get_new_depth_vel(reshape(model,(/nv(1),nv(2)/)),tau_pars& &,vel_pars) !################ !WRITE OUT MODEL# !################ call srite("out",vel_pars%velz,size(vel_pars%velz)*4) end if end if end program