program MatchRay c----------------------------------------------- c c MatchRay c c----------------------------------------------------------------- c Purpose: Compute ray parameters along rays in a model c-------------------------------------------------------------------------- c Input file: c c------------------------------------------------------------------------ c Output file: Gives ray paths that can be plotted with Graph. c c------------------------------------------------------------------------ c% implicit none include '../Src/rayparz.inc' logical ray_in_tau,found,found_multi,end_ray,ray_far_ref integer esize,esize4,it,idelt,iref,one,count_match integer fetch,getch,auxpar,sreed,sseek,srite,auxin,tempi integer maxt,ndelt,nsou,nx_ref,n_ref,isou,ix_ref,n_mid,n_h integer nby_ray,nby_seek,nby_sou,nby_mod integer idelt_ref,idelt_inc,idelt_ref_min,idelt_inc_min integer isou_ref,isou_inc,imid,ih integer nt_inc,nt_ref,zmap,n_zmap integer want_multi,cross_time,cross_only,cross_dir,bad integer it_min,ix_ref_min,ix_ref_max,dipfd real xray1,xray2,vray1,vray2,xref1,xref2,vref1,vref2,vref0 real wxref,xref,vref,pxref,pvref,velref,sigref real xa1,xa2,xb1,xb2,pa1,pa2,pb1,pb2,va1,vb1 real dx_ref,idx_ref,tmatch_ref,tmatch_inc,tmatch_tot real mid_0,h_0,d_mid,d_h,sou_inc,sou_ref real wt,w_inc,w_ref real sin_cross,cos_cross,tandip_cross real tmin_ref,vmin_ref,vmax_ref real delta,delt,delmin real x0_ref,dip_ang real xsou,zsou,tausou,dx_sou,dz_sou,dtau_sou double precision depsi,sepsi,det,a,b,a_int,b_int double precision Dxa,Dva,Dpa,Dxb,Dvb,Dpb double precision Dxab,Dvab,Dpab logical crossed(mangle,msou) real rtraj(2,mtim,mangle),rpar(2,mtim,mangle) real rmod(2,mtim,mangle) real rtraj_inc(2,mtim,2),rpar_inc(2,mtim,2) real rmod_inc(2,mtim,2) real rtraj_ref(2,mtim,2),rpar_ref(2,mtim,2) real rmod_ref(2,mtim,2) real rtraj_zmap_inc(2,mtim,2),rtraj_zmap_ref(2,mtim,2) real rtraj_match(2,2*mtim),rpar_match(2,2) real rmod_match(2,2*mtim) real refl(2,mnx_ref) real cosrefl(mnx_ref),sinrefl(mnx_ref),tandip(mnx_ref) real tcross(mangle,msou) real xvcross(2,mangle,msou) real pcross_inc(2,mangle,msou) real dipcross(mangle,msou) real pcross_ref(2,mangle,msou) real velcross_inc(mangle,msou) real velcross_ref(mangle,msou) real sigcross_inc(mangle,msou) real sigcross_ref(mangle,msou) real location(3) call initpar() call doc( * '/home/kana/sep/biondo/Code/Focus/Raytrace/MatchRay/MatchRay.f') c set constants depsi= 1.11D-15 sepsi= 1.0d-4 c fetch dimensions and esize. if (fetch('n1','d',maxt).eq.0) * call erexit('n1 must be supplied') if (fetch('n2','d',ndelt).eq.0) * call erexit('n2 must be supplied') if (fetch('d1','f',delt).eq.0) * call erexit('d1 must be supplied') if (fetch('d2','f',delta).eq.0) * call erexit('d2 must be supplied') if (fetch('o2','f',delmin).eq.0) * call erexit('o2 must be supplied') if (fetch('nsou','d',nsou).eq.0) then if(fetch('n3','d',nsou).eq.0) * call erexit('nsou must be supplied') else ndelt=ndelt/nsou end if if (fetch('esize','d',esize).eq.0) * call erexit('esize must be supplied') c get parametes from parameter file if (getch('tmin_ref','f',tmin_ref).le.0) tmin_ref=0. call putch('tmin_ref','f',tmin_ref) it_min=int(tmin_ref/delt) write(0,*)'it_min=',it_min if (getch('zmap','d',zmap).le.0) zmap=0 call putch('zmap','d',zmap) if(zmap .eq. 1) then if (auxpar('n1','d',n_zmap,'Rayzmap').eq.0) * call erexit('n1 must be supplied in Rayzmap') end if if (getch('want_multi','d',want_multi).le.0) want_multi=0 call putch('want_multi','d',want_multi) if (getch('cross_time','d',cross_time).le.0) cross_time=0 call putch('cross_time','d',cross_time) if (getch('cross_only','d',cross_only).le.0) cross_only=0 call putch('cross_only','d',cross_only) if (getch('cross_dir','d',cross_dir).le.0) cross_dir=1 call putch('cross_dir','d',cross_dir) if(cross_only .ne. 1) then c { if (fetch('xsou','f',xsou).eq.0) * call erexit('xsou must be supplied') if (fetch('dx_sou','f',dx_sou).eq.0) * call erexit('dx_sou must be supplied') if (fetch('dz_sou','f',dz_sou).eq.0) then if (fetch('dtau_sou','f',dtau_sou).eq.0) then if(fetch('output_tau','d',tempi).eq.0) then call erexit('Either dz_sou or dtau_sou or output_tau must be supplied') else if(tempi.eq.1) ray_in_tau=.true. end if else ray_in_tau=.true. end if end if c if (fetch('tausou','f',tausou).eq.0) c * call erexit('tausou must be supplied') c else c if (fetch('zsou','f',zsou).eq.0) c * call erexit('zsou must be supplied') c ray_in_tau=.false. c end if write(0,*)'ray_in_tau=',ray_in_tau if (getch('n_mid','d',n_mid).eq.0) n_mid=1 if (getch('mid_0','f',mid_0).eq.0) mid_0=xsou if (getch('d_mid','f',d_mid).eq.0) d_mid=dx_sou if (getch('n_h','d',n_h).eq.0) n_h=1 if (getch('h_0','f',h_0).eq.0) h_0=0. if (getch('d_h','f',d_h).eq.0) d_h=dx_sou if ((mid_0+(n_mid-1)*d_mid + h_0+(n_h-1)*d_h) * .gt. (xsou+(nsou-1)*dx_sou)) * call erexit('out on the far side') if ((mid_0 - h_0+(n_h-1)*d_h) .lt. xsou) * call erexit('out on the near side') c } end if C get information on reflector if (auxpar('n1','d',nx_ref,'Refl').eq.0) * call erexit('n1 must be supplied in Refl') if (auxpar('n2','d',n_ref,'Refl').eq.0) * call erexit('n2 must be supplied in Refl') C if (auxpar('o1','f',x0_ref,'Refl').eq.0) C * call erexit('o1 must be supplied in Refl') C if (auxpar('d1','f',dx_ref,'Refl').eq.0) C * call erexit('d1 must be supplied in Refl') c check parameters if (ndelt .gt. mangle) call erexit('ndelt too large') if (maxt .gt. mtim) call erexit('maxt too large') if (nsou .gt. msou) call erexit('nsou too large') if (nx_ref .gt. mnx_ref) call erexit('nx_ref too large') if(cross_time .eq. 1) then if(nx_ref .ne. nsou*ndelt) then write(0,*) "nx_ref=",nx_ref write(0,*) "nsou=",nsou write(0,*) "ndelt=",ndelt call erexit('nx_ref != nsou*ndelt') end if end if C if (mod(mid_0-xsou,dx_sou) .ne. 0) C * call erexit('mid_0-xsou not multiple of dx_sou') C if (mod(d_mid,dx_sou) .ne. 0) C * call erexit('d_mid not multiple of dx_sou') C if (mod(d_h,dx_sou) .ne. 0) C * call erexit('d_h not multiple of dx_sou') esize4=4 one=1 call putch('esize','d',esize) call putch('n1','d',2*maxt) call putch('o1','f',0.) call putch('d1','f',delt) if(want_multi .ne. 1) then call putch('n2','d',n_h) call putch('n3','d',n_mid) call putch('n4','d',n_ref) call putch('o2','f',h_0) call putch('o3','f',mid_0) call putch('o4','f',1.) call putch('d2','f',d_h) call putch('d3','f',d_mid) call putch('d4','f',1.) end if call putch('n_mid','d',n_mid) call putch('mid_0','f',mid_0) call putch('d_mid','f',d_mid) call putch('n_h','d',n_h) call putch('h_0','f',h_0) call putch('d_h','f',d_h) call putch('n_ref','d',n_ref) call putch('ref_0','f',1.) call putch('d_ref','f',1.) if(cross_only .ne. 1) then c { call auxputch('esize','d',esize,'Raypar_match') call auxputch('esize','d',esize4,'Time_match') call auxputch('esize','d',esize4,'Time_ref') call auxputch('esize','d',esize,'Raymod_match') call auxputch('n1','d',2,'Raypar_match') call auxputch('n1','d',2*maxt,'Raymod_match') call auxputch('o1','f',0.,'Raypar_match') call auxputch('o1','f',0.,'Raymod_match') call auxputch('d1','f',delt,'Raypar_match') call auxputch('d1','f',delt,'Raymod_match') if(want_multi .ne. 1) then call auxputch('n1','d',n_h,'Time_match') call auxputch('n1','d',n_h,'Time_ref') call auxputch('n2','d',n_h,'Raymod_match') call auxputch('n2','d',n_h,'Raypar_match') call auxputch('n2','d',n_mid,'Time_match') call auxputch('n2','d',n_mid,'Time_ref') call auxputch('n2','d',n_mid,'Pos_match') call auxputch('n3','d',n_mid,'Raymod_match') call auxputch('n3','d',n_mid,'Raypar_match') call auxputch('n3','d',n_ref,'Time_match') call auxputch('n3','d',n_ref,'Time_ref') call auxputch('n4','d',n_ref,'Raypar_match') call auxputch('n4','d',n_ref,'Raymod_match') call auxputch('o1','f',h_0,'Time_match') call auxputch('o1','f',h_0,'Time_ref') call auxputch('o2','f',h_0,'Raymod_match') call auxputch('o2','f',h_0,'Raypar_match') call auxputch('o2','f',mid_0,'Time_match') call auxputch('o2','f',mid_0,'Time_ref') call auxputch('o3','f',mid_0,'Raypar_match') call auxputch('o3','f',mid_0,'Raymod_match') call auxputch('o4','f',one,'Raypar_match') call auxputch('o4','f',one,'Raymod_match') call auxputch('o3','f',one,'Time_match') call auxputch('o3','f',one,'Time_ref') call auxputch('d1','f',d_h,'Time_match') call auxputch('d1','f',d_h,'Time_ref') call auxputch('d2','f',d_h,'Raypar_match') call auxputch('d2','f',d_h,'Raymod_match') call auxputch('d2','f',d_mid,'Time_match') call auxputch('d2','f',d_mid,'Time_ref') call auxputch('d3','f',d_mid,'Raypar_match') call auxputch('d3','f',d_mid,'Raymod_match') call auxputch('d3','f',1.,'Time_match') call auxputch('d3','f',1.,'Time_ref') call auxputch('d4','f',1.,'Raypar_match') call auxputch('d4','f',1.,'Raymod_match') else call auxputch('n1','d',3,'location') end if c } else c { call auxputch('esize','d',esize4,'Time_cross') call auxputch('esize','d',esize,'Xv_cross') call auxputch('esize','d',esize,'P_cross') call auxputch('n1','d',nsou,'Time_cross') call auxputch('n2','d',n_ref,'Time_cross') call auxputch('n3','d',one,'Time_cross') call auxputch('d1','f',1.,'Time_cross') call auxputch('d2','f',1.,'Time_cross') call auxputch('d3','f',1.,'Time_cross') call auxputch('o1','f',0.,'Time_cross') call auxputch('o2','f',0.,'Time_cross') call auxputch('o3','f',0.,'Time_cross') call auxputch('n1','d',nsou,'Xv_cross') call auxputch('n2','d',n_ref,'Xv_cross') call auxputch('n3','d',one,'Xv_cross') call auxputch('d1','f',1.,'Xv_cross') call auxputch('d2','f',1.,'Xv_cross') call auxputch('d3','f',1.,'Xv_cross') call auxputch('o1','f',0.,'Xv_cross') call auxputch('o2','f',0.,'Xv_cross') call auxputch('o3','f',0.,'Xv_cross') call auxputch('n1','d',nsou,'P_cross') call auxputch('n2','d',n_ref,'P_cross') call auxputch('n3','d',one,'P_cross') call auxputch('d1','f',1.,'P_cross') call auxputch('d2','f',1.,'P_cross') call auxputch('d3','f',1.,'P_cross') call auxputch('o1','f',0.,'P_cross') call auxputch('o2','f',0.,'P_cross') call auxputch('o3','f',0.,'P_cross') call auxputch('n1','d',nsou,'Dip_cross') call auxputch('n2','d',n_ref,'Dip_cross') call auxputch('n3','d',one,'Dip_cross') call auxputch('d1','f',1.,'Dip_cross') call auxputch('d2','f',1.,'Dip_cross') call auxputch('d3','f',1.,'Dip_cross') call auxputch('o1','f',0.,'Dip_cross') call auxputch('o2','f',0.,'Dip_cross') call auxputch('o3','f',0.,'Dip_cross') call auxputch('esize','d',8,'Pos_match') call auxputch('n1','d',n_h,'Pos_match') call auxputch('n3','d',n_ref,'Pos_match') call auxputch('o1','f',h_0,'Pos_match') call auxputch('o3','f',mid_0,'Pos_match') call auxputch('d1','f',d_h,'Pos_match') call auxputch('d2','f',d_mid,'Pos_match') call auxputch('d3','f',1.,'Pos_match') end if c } c close history file bad=0 count_match=0 C call hclose() C call auxclose('Raypar_match') c precompute information on end point of each ray C idx_ref=1./dx_ref nby_ray=maxt*esize nby_sou=nby_ray*ndelt dipfd=auxin('Dip_refl') do iref=1,n_ref c { location(3)=iref c read reflector information if (sreed('Refl',refl,nx_ref*8) .ne. nx_ref*8) then call erexit("Problems reading Refl file") end if write(0,*)'Matching rays for reflector #',iref vmin_ref=100000000. vmax_ref=-100000000. do ix_ref=1,nx_ref c { if(vmin_ref .gt. refl(2,ix_ref)) then vmin_ref=refl(2,ix_ref) end if if(vmax_ref .lt. refl(2,ix_ref)) then vmax_ref=refl(2,ix_ref) end if c } end do write(0,*)'vmin_ref=',vmin_ref write(0,*)'vmax_ref=',vmax_ref if(dipfd .gt. 0) then c { write(0,*)'read Dip_refl file' write(0,*)'dipfd ',dipfd if (sreed('Dip_refl',tandip,nx_ref*4) .ne. nx_ref*4) then call erexit("Problems reading Dip_refl file") end if c } else c { write(0,*)'NOT read Dip_refl file' c comute cos and sin of 2*dip angle do ix_ref=2,nx_ref-1 c { tandip(ix_ref)= * (refl(2,ix_ref+1)-refl(2,ix_ref-1))/ * (refl(1,ix_ref+1)-refl(1,ix_ref-1)) end do c } tandip(1) = * (refl(2,2)-refl(2,1))/ * (refl(1,2)-refl(1,1)) tandip(nx_ref) = * (refl(2,nx_ref)-refl(2,nx_ref-1))/ * (refl(1,nx_ref)-refl(1,nx_ref-1)) c } end if do ix_ref=1,nx_ref c { dip_ang=2*atan(tandip(ix_ref)) sinrefl(ix_ref)=sin(dip_ang) cosrefl(ix_ref)=cos(dip_ang) c } end do call snap("Tandip.h",nx_ref,1,tandip) do isou=1,nsou c { do idelt=1,ndelt c { if(iref .ne. 1) then nby_seek=nby_sou*(isou-1) + nby_ray*(idelt-1) if (sseek('in',nby_seek,0).ne.nby_seek) then call erexit("Problems seeking Input file") end if if (sseek('Raypar',nby_seek,0).ne.nby_seek) then call erexit("Problems seeking Raypar file") end if if (sseek('Raymod',nby_seek,0).ne.nby_seek) then call erexit("Problems seeking Raymod file") end if end if if (sreed('in',rtraj(1,1,idelt),nby_ray).ne.nby_ray) then call erexit("Problems reading Input file") end if if (sreed('Raypar',rpar(1,1,idelt),nby_ray).ne.nby_ray) then call erexit("Problems reading Raypar file") end if if (sreed('Raymod',rmod(1,1,idelt),nby_ray).ne.nby_ray) then call erexit("Problems reading Raymod file") end if c if(isou.eq.1 .and. iref.ne.1) then c if (sseek('in',0,0).ne.0) then c call erexit("Problems seeking Input file") c end if c if (sseek('Raypar',0,0).ne.0) then c call erexit("Problems seeking Raypar file") c end if c if (sseek('Raymod',0,0).ne.0) then c call erexit("Problems seeking Raymod file") c end if c end if c c if (sreed('in',rtraj(1,1,1),nby_sou).ne.nby_sou) then c call erexit("Problems reading Input file") c end if c c if (sreed('Raypar',rpar(1,1,1),nby_sou).ne.nby_sou) then c call erexit("Problems reading Raypar file") c end if c cc if (sreed('Raymod',rmod(1,1,1),nby_sou).ne.nby_sou) then c call erexit("Problems reading Raymod file") c end if c c do idelt=1,ndelt end_ray=.false. crossed(idelt,isou) = .false. it=it_min do while ((.not.crossed(idelt,isou)).and. * (it.lt.maxt-1).and. * (.not. end_ray)) c { it=it+1 xray1=rtraj(1,it ,idelt) xray2=rtraj(1,it+1,idelt) vray1=rtraj(2,it ,idelt) vray2=rtraj(2,it+1,idelt) if((vray1 .eq. vray2).and.(xray1 .eq. xray1)) then end_ray=.true. end if if(cross_time .eq. 1) then vray1=(it-1.)*delt vray2=(it )*delt ix_ref_min=idelt+(isou-1)*ndelt-1 ix_ref_max=ix_ref_min+1 else ix_ref_min=0 ix_ref_max=nx_ref-1 end if Dxa=xray2-xray1 Dva=vray2-vray1 if (cross_dir .eq. 1) then if((vray1 .lt. vmin_ref).and.(vray2 .lt. vmin_ref)) then ray_far_ref=.true. else ray_far_ref=.false. end if else if((vray1 .gt. vmax_ref).and.(vray2 .gt. vmax_ref)) then ray_far_ref=.true. else ray_far_ref=.false. end if end if ix_ref=ix_ref_min do while ((.not.crossed(idelt,isou)).and. * (ix_ref.lt.ix_ref_max).and. * (.not. end_ray).and. * (.not. ray_far_ref)) c { ix_ref=ix_ref+1 c trick to keep the same code for both cases (more expensive) if(cross_time .eq. 1) then xref1=xray1-sign(sepsi,Dxa) xref2=xray2+sign(sepsi,Dxa) vref1=refl(2,ix_ref) vref2=vref1 else xref1=refl(1,ix_ref ) xref2=refl(1,ix_ref+1) vref1=refl(2,ix_ref ) vref2=refl(2,ix_ref+1) end if Dxb=xref2-xref1 Dvb=vref2-vref1 Dxab=xref1-xray1 Dvab=vref1-vray1 det=Dva*Dxb-Dxa*Dvb if(dabs(det) .le. depsi) then c { if((Dxa .eq. Dxb) .and. (Dva .eq. Dvb) .and. * (Dxab .eq. Dvab)) then a=0.d0 b=0.d0 call seperr('crossing ray is coincident with reflector') else write(0,*)'POTENTIAL ERROR' write(0,*)'ray is parallel to reflctor' write(0,*)'det=',det write(0,*)'it=',it write(0,*)'ix_ref=',ix_ref write(0,*)'isou=',isou write(0,*)'idelt=',idelt write(0,*)'xray1=',xray1 write(0,*)'xray2=',xray2 write(0,*)'xref1=',xref1 write(0,*)'xref2=',xref2 write(0,*)'vray1=',vray1 write(0,*)'vray2=',vray2 write(0,*)'vref1=',vref1 write(0,*)'vref2=',vref2 write(0,*)'Dxab=',Dxab write(0,*)'Dvab=',Dvab C call seperr('ray is parallel to reflctor') end if c } else c { a= (Dvab*Dxb - Dxab*Dvb)/det b=-(-Dvab*Dxa + Dxab*Dva)/det if((a .ge. -sepsi) .and. (a .le. 1.+sepsi) .and. * (b .ge. -sepsi) .and. (b .le. 1.+sepsi)) then c { crossed(idelt,isou) = .true. xvcross(1,idelt,isou)=xray1 + (a*Dxa) if(cross_time .eq. 1) then xvcross(2,idelt,isou)=rtraj(2,it ,idelt)*(1.-a) + * rtraj(2,it+1,idelt)*( a) else xvcross(2,idelt,isou)=vray1 + (a*Dva) end if tcross(idelt,isou)=((it-1.)+a)*delt C write(0,*)'xvcross(1)=',xvcross(1,idelt,isou) C write(0,*)'xvcross(2)=',xvcross(2,idelt,isou) C write(0,*)'tcross(2)=',tcross(idelt,isou) pcross_inc(1,idelt,isou) = * rpar(1,it ,idelt)*(1.-a) + * rpar(1,it+1,idelt)*( a) pcross_inc(2,idelt,isou) = * rpar(2,it ,idelt)*(1.-a) + * rpar(2,it+1,idelt)*( a) velcross_inc(idelt,isou) = * rmod(1,it ,idelt)*(1.-a) + * rmod(1,it+1,idelt)*( a) sigcross_inc(idelt,isou) = * rmod(2,it ,idelt)*(1.-a) + * rmod(2,it+1,idelt)*( a) dipcross(idelt,isou) = * -tan(asin((pcross_inc (1,idelt,isou) * * velcross_inc( idelt,isou)))) if(cross_only .ne. 1) then c { c note that they have the opposite sign that they should because c all rays are going down if(.not. ray_in_tau) then c { sin_cross= * sinrefl(ix_ref )*(1.-b) + * sinrefl(ix_ref+1)*( b) cos_cross= * cosrefl(ix_ref )*(1.-b) + * cosrefl(ix_ref+1)*( b) pcross_ref(1,idelt,isou) = * pcross_inc(1,idelt,isou)*cos_cross + * pcross_inc(2,idelt,isou)*sin_cross pcross_ref(2,idelt,isou) = * pcross_inc(1,idelt,isou)*sin_cross - * pcross_inc(2,idelt,isou)*cos_cross c } else c { if(dipfd .le. 0) then c { tandip_cross= * tandip(ix_ref )*(1.-b) + * tandip(ix_ref+1)*( b) C ATTENTION VALID ONLY FOR FLAT REFLECTORS dip_ang=2*atan(.5*velcross_inc(idelt,isou) * * (tandip_cross-sigcross_inc(idelt,isou))) sin_cross=sin(dip_ang) cos_cross=cos(dip_ang) c } else c { c in this case the angles are already in depth and not in tau sin_cross= * sinrefl(ix_ref )*(1.-b) + * sinrefl(ix_ref+1)*( b) cos_cross= * cosrefl(ix_ref )*(1.-b) + * cosrefl(ix_ref+1)*( b) c } end if pcross_ref(1,idelt,isou) = * pcross_inc(1,idelt,isou)* * (cos_cross - * sin_cross* * .5*sigcross_inc(idelt,isou)*velcross_inc(idelt,isou)) + * pcross_inc(2,idelt,isou)* * ( * cos_cross* * 2.*sigcross_inc(idelt,isou) + * sin_cross* * (2./velcross_inc(idelt,isou) - * .5*velcross_inc(idelt,isou)*sigcross_inc(idelt,isou)**2) * ) pcross_ref(2,idelt,isou) = * +pcross_inc(1,idelt,isou)* * (sin_cross* * .5*velcross_inc(idelt,isou)) + * pcross_inc(2,idelt,isou)* * ( * -cos_cross + * sin_cross* * .5*sigcross_inc(idelt,isou)*velcross_inc(idelt,isou) * ) c } end if c } end if c } end if c } end if c } end do c } end do if (.not. crossed(idelt,isou)) then xvcross(1,idelt,isou)=rtraj(1,1,idelt) xvcross(2,idelt,isou)=rtraj(2,1,idelt) tcross(idelt,isou)=-1. if(cross_time .eq. 1) then write(0,*)'isou=',isou write(0,*)'idelt=',idelt call seperr('Ray never crossed time surface') end if end if c } end do c } end do write(0,*) "finished reflector" do isou=1,nsou idelt= srite("xvcross.H",xvcross(1,1,isou),ndelt*8) idelt= srite("tcross.H",tcross(1,isou),ndelt*4) end do if(cross_only .eq. 1) then c { do isou=1,nsou c { if (srite('Time_cross',tcross(1,isou),ndelt*4) * .ne. ndelt*4) then call erexit("Problems writing Time_cross file") end if if (srite('Xv_cross',xvcross(1,1,isou),ndelt*8) * .ne. ndelt*8) then call erexit("Problems writing Xv_cross file") end if if (srite('P_cross',pcross_inc(1,1,isou),ndelt*8) * .ne. ndelt*8) then call erexit("Problems writing P_cross file") end if if (srite('Dip_cross',dipcross(1,isou),ndelt*4) * .ne. ndelt*4) then call erexit("Problems writing Dip_cross file") end if c } end do c } else c { write(0,*)'BEGINNING matching rays' C call snap("Xcross",mangle,nsou,xcross) C call snap("Tcross",mangle,nsou,tcross) C call snap("Vcross",mangle,nsou,vcross) C call snap("Pxcross_inc",mangle,nsou,pxcross_inc) C call snap("Pvcross_inc",mangle,nsou,pvcross_inc) C call snap("Pxcross_ref",mangle,nsou,pxcross_ref) C call snap("Pvcross_ref",mangle,nsou,pvcross_ref) found_multi=.false. do imid=1,n_mid c { location(2)=imid do ih=1,n_h c { location(1)=ih sou_inc=(mid_0 + (imid-1)*d_mid) - (h_0 + (ih-1)*d_h) - xsou sou_ref=(mid_0 + (imid-1)*d_mid) + (h_0 + (ih-1)*d_h) - xsou C if(abs(mod(sou_inc,dx_sou)) .gt. dx_sou*.001) then C write(0,*)'mod(sou_inc,dx_sou)=',mod(sou_inc,dx_sou) C write(0,*)'sou_inc=',sou_inc C call seperr('sou_inc not on grid') C else C write(0,*)'sou_inc=',sou_inc isou_inc= nint(sou_inc/dx_sou) + 1 C end if C if(abs(mod(sou_ref,dx_sou)) .gt. dx_sou*.001) then C write(0,*)'mod(sou_ref,dx_sou)=',mod(sou_ref,dx_sou) C write(0,*)'sou_ref=',sou_ref C call seperr('sou_ref not on grid') C else C write(0,*)'sou_ref=',sou_ref isou_ref= nint(sou_ref/dx_sou) + 1 C end if C write(0,*)'imid=',imid C write(0,*)'ih=',ih C write(0,*)'isou_inc=',isou_inc C write(0,*)'isou_ref=',isou_ref idelt_inc=0 found=.false. C do while ((idelt_inc .lt. ndelt-1).and.(.not. found)) do while ((idelt_inc .lt. ndelt-1)) c { idelt_inc=idelt_inc+1 if(crossed(idelt_inc ,isou_inc) .and. * crossed(idelt_inc+1,isou_inc)) then c { idelt_ref=0 C do while ((idelt_ref .lt. ndelt-1).and.(.not. found)) do while ((idelt_ref .lt. ndelt-1)) c { idelt_ref=idelt_ref+1 if(crossed(idelt_ref ,isou_ref) .and. * crossed(idelt_ref+1,isou_ref)) then c { xa1=xvcross(1,idelt_inc ,isou_inc) xa2=xvcross(1,idelt_inc+1,isou_inc) xb1=xvcross(1,idelt_ref ,isou_ref) xb2=xvcross(1,idelt_ref+1,isou_ref) pa1=-pcross_ref(1,idelt_inc ,isou_inc) pa2=-pcross_ref(1,idelt_inc+1,isou_inc) pb1=pcross_inc(1,idelt_ref ,isou_ref) pb2=pcross_inc(1,idelt_ref+1,isou_ref) Dxa=xa2-xa1 Dpa=pa2-pa1 Dxb=xb2-xb1 Dpb=pb2-pb1 Dxab=xb1-xa1 Dpab=pb1-pa1 det=Dpa*Dxb-Dxa*Dpb if(dabs(det) .le. depsi) then c { if((Dxa .eq. Dxb) .and. (Dpa .eq. Dpb) .and. * (Dxab .eq. Dpab)) then a=0.d0 b=0.d0 found=.true. call seperr('matching rays are same') else det=depsi write(0,*)'ERROR' write(0,*)'idelt_inc=',idelt_inc write(0,*)'idelt_ref=',idelt_ref write(0,*)'xa1=',xa1 write(0,*)'xa2=',xa2 write(0,*)'xb1=',xb1 write(0,*)'xb2=',xb2 write(0,*)'Dxab=',Dxab write(0,*)'Dpab=',Dpab call seperr('matching rays are parrallel') end if c } else c { a= (Dpab*Dxb - Dxab*Dpb)/det b=-(-Dpab*Dxa + Dxab*Dpa)/det if((a .ge. -sepsi) .and. (a .le. 1.+sepsi) .and. * (b .ge. -sepsi) .and. (b .le. 1.+sepsi)) then c { C write(0,*)'GOOD det=',det C write(0,*)'GOOD a=',a C write(0,*)'GOOD b=',b C write(0,*)'idelt_inc=',idelt_inc C write(0,*)'idelt_ref=',idelt_ref C write(0,*)'xa1=',xa1 C write(0,*)'xa2=',xa2 C write(0,*)'xb1=',xb1 C write(0,*)'xb2=',xb2 C write(0,*)'pa1=',pa1 C write(0,*)'pa2=',pa2 C write(0,*)'pb1=',pb1 C write(0,*)'pb2=',pb2 C write(0,*) "ANGLE CROSSS?",pcross_ref(1,idelt_inc ,isou_inc) C * ,pcross_ref(1,idelt_inc +1,isou_inc) C * ,pcross_inc(1,idelt_ref ,isou_ref) C * ,pcross_inc(1,idelt_ref +1,isou_ref) C if(found .and. (want_multi .ne. 1)) then c { found_multi=.true. c } else c { if(found) found_multi=.true. C write(0,*) "xa1 xa2 dxa",xa1,xa2,dxa C write(0,*) "pa1 pa2 dpa",pa1,pa2,dpa C write(0,*) "xb1 xb2 dxb",xa2,xa2,dxa C write(0,*) "pb1 pb2 dpb",pb1,pb2,dpb C write(0,*) "Dxab Dpab det",Dxab,Dpab,det found=.true. count_match=count_match+1 xref=xa1+ (a*Dxa) vref= * xvcross(2,idelt_inc ,isou_inc)*(1.-a) + * xvcross(2,idelt_inc+1,isou_inc)*(a ) pxref= * pcross_inc(1,idelt_inc ,isou_inc)*(1.-a) + * pcross_inc(1,idelt_inc+1,isou_inc)*(a ) pvref= * pcross_inc(2,idelt_inc ,isou_inc)*(1.-a) + * pcross_inc(2,idelt_inc+1,isou_inc)*(a ) velref= * velcross_inc(idelt_inc ,isou_inc)*(1.-a) + * velcross_inc(idelt_inc+1,isou_inc)*(a ) sigref= * sigcross_inc(idelt_inc ,isou_inc)*(1.-a) + * sigcross_inc(idelt_inc+1,isou_inc)*(a ) tmatch_inc= * tcross(idelt_inc ,isou_inc)*(1.-a) + * tcross(idelt_inc+1,isou_inc)*(a ) tmatch_ref= * tcross(idelt_ref ,isou_ref)*(1.-b) + * tcross(idelt_ref+1,isou_ref)*(b ) tmatch_tot=tmatch_inc+tmatch_ref C write(0,*)'xref=',xref C write(0,*)'vref=',vref C write(0,*)'pxref=',pxref C write(0,*)'pvref=',pvref C write(0,*)'velref=',velref C write(0,*)'tmatch_tot=',tmatch_tot C write(0,*)'tcross(idelt_inc)=', C * tcross(idelt_inc,isou_inc) C write(0,*)'tcross(idelt_inc+1)=', C * tcross(idelt_inc+1,isou_inc) C write(0,*)'tcross(idelt_ref)=', C * tcross(idelt_ref,isou_ref) C write(0,*)'tcross(idelt_ref+1)=', C * tcross(idelt_ref+1,isou_ref) C C write(0,*)'xcross(idelt_inc)=', C * xcross(idelt_inc,isou_inc) C write(0,*)'xcross(idelt_inc+1)=', C * xcross(idelt_inc+1,isou_inc) C write(0,*)'xcross(idelt_ref)=', C * xcross(idelt_ref,isou_ref) C write(0,*)'xcross(idelt_ref+1)=', C * xcross(idelt_ref+1,isou_ref) C C write(0,*)'pxcross_ref(idelt_inc)=', C * pxcross_ref(idelt_inc,isou_inc) C write(0,*)'pxcross_ref(idelt_inc+1)=', C * pxcross_ref(idelt_inc+1,isou_inc) C write(0,*)'pxcross_inc(idelt_inc)=', C * pxcross_inc(idelt_inc,isou_inc) C write(0,*)'pxcross_inc(idelt_inc+1)=', C * pxcross_inc(idelt_inc+1,isou_inc) nby_seek=nby_sou*(isou_inc-1) + nby_ray*(idelt_inc-1) if (sseek('in',nby_seek,0).ne.nby_seek) then call erexit("Problems seeking Input file") end if if (sreed('in',rtraj_inc(1,1,1),nby_ray).ne.nby_ray) then call erexit("Problems reading Input file") end if if (sreed('in',rtraj_inc(1,1,2),nby_ray).ne.nby_ray) then call erexit("Problems reading Input file") end if if (sseek('Raypar',nby_seek,0).ne.nby_seek) then call erexit("Problems seeking Raypar file") end if if (sreed('Raypar',rpar_inc(1,1,1),nby_ray).ne.nby_ray) then call erexit("Problems reading Raypar file") end if if (sreed('Raypar',rpar_inc(1,1,2),nby_ray).ne.nby_ray) then call erexit("Problems reading Raypar file") end if if (sseek('Raymod',nby_seek,0).ne.nby_seek) then call erexit("Problems seeking Raymod file") end if if (sreed('Raymod',rmod_inc(1,1,1),nby_ray).ne.nby_ray) then call erexit("Problems reading Raymod file") end if if (sreed('Raymod',rmod_inc(1,1,2),nby_ray).ne.nby_ray) then call erexit("Problems reading Raymod file") end if if(zmap .eq. 1) then if (sseek('Rayzmap',nby_seek,0).ne.nby_seek) then call erexit("Problems seeking Rayzmap file") end if if (sreed('Rayzmap',rtraj_zmap_inc(1,1,1),nby_ray) * .ne.nby_ray) then call erexit("Problems reading Rayzmap file") end if if (sreed('Rayzmap',rtraj_zmap_inc(1,1,2),nby_ray) * .ne.nby_ray) then call erexit("Problems reading Rayzmap file") end if end if nby_seek=nby_sou*(isou_ref-1) + nby_ray*(idelt_ref-1) if (sseek('in',nby_seek,0).ne.nby_seek) then call erexit("Problems seeking Input file") end if if (sreed('in',rtraj_ref(1,1,1),nby_ray).ne.nby_ray) then call erexit("Problems reading Input file") end if if (sreed('in',rtraj_ref(1,1,2),nby_ray).ne.nby_ray) then call erexit("Problems reading Input file") end if if (sseek('Raypar',nby_seek,0).ne.nby_seek) then call erexit("Problems seeking Raypar file") end if if (sreed('Raypar',rpar_ref(1,1,1),nby_ray).ne.nby_ray) then call erexit("Problems reading Raypar file") end if if (sreed('Raypar',rpar_ref(1,1,2),nby_ray).ne.nby_ray) then call erexit("Problems reading Raypar file") end if if (sseek('Raymod',nby_seek,0).ne.nby_seek) then call erexit("Problems seeking Raymod file") end if if (sreed('Raymod',rmod_ref(1,1,1),nby_ray).ne.nby_ray) then call erexit("Problems reading Raymod file") end if if (sreed('Raymod',rmod_ref(1,1,2),nby_ray).ne.nby_ray) then call erexit("Problems reading Raymod file") end if if(zmap .eq. 1) then if (sseek('Rayzmap',nby_seek,0).ne.nby_seek) then call erexit("Problems seeking Rayzmap file") end if if (sreed('Rayzmap',rtraj_zmap_ref(1,1,1),nby_ray) * .ne.nby_ray) then call erexit("Problems reading Rayzmap file") end if if (sreed('Rayzmap',rtraj_zmap_ref(1,1,2),nby_ray) * .ne.nby_ray) then call erexit("Problems reading Rayzmap file") end if end if nt_inc=(tmatch_inc)/delt+1 nt_ref=(tmatch_ref)/delt+1 if((nt_inc+nt_ref +1) .gt. 2*maxt) then call seperr('maxt too small') end if wt=tmatch_inc/delt-(nt_inc-1.) C wt=1. xa1= rtraj_inc(1,nt_inc ,1)*(1.-wt) + * rtraj_inc(1,nt_inc+1,1)*( wt) xb1= rtraj_inc(1,nt_inc ,2)*(1.-wt) + * rtraj_inc(1,nt_inc+1,2)*( wt) w_inc=(xref-xa1)/(xb1-xa1) if(zmap .eq. 1) then va1= rtraj_zmap_inc(2,nt_inc ,1)*(1.-wt) + * rtraj_zmap_inc(2,nt_inc+1,1)*( wt) vb1= rtraj_zmap_inc(2,nt_inc ,2)*(1.-wt) + * rtraj_zmap_inc(2,nt_inc+1,2)*( wt) vref=va1+w_inc*(vb1-va1) end if C write(0,*)'tmatch_inc',tmatch_inc C write(0,*)'wt',wt C write(0,*)'xa1',xa1 C write(0,*)'xb1',xb1 C write(0,*)'va1',pa1 C write(0,*)'vb1',pb1 C write(0,*)'xref',xref C write(0,*)'idelt_inc',idelt_inc C write(0,*)'idelt_ref',idelt_ref C write(0,*)'VREF',vref C wt=tmatch_ref/delt-(nt_ref-1.) xa1= rtraj_ref(1,nt_ref ,1)*(1.-wt) + * rtraj_ref(1,nt_ref+1,1)*( wt) xb1= rtraj_ref(1,nt_ref ,2)*(1.-wt) + * rtraj_ref(1,nt_ref+1,2)*( wt) w_ref=(xref-xa1)/(xb1-xa1) if(zmap .eq. 1) then va1= rtraj_zmap_ref(2,nt_ref ,1)*(1.-wt) + * rtraj_zmap_ref(2,nt_ref+1,1)*( wt) vb1= rtraj_zmap_ref(2,nt_ref ,2)*(1.-wt) + * rtraj_zmap_ref(2,nt_ref+1,2)*( wt) vref=va1+w_ref*(vb1-va1) end if C write(0,*)'tmatch_ref',tmatch_ref C write(0,*)'wt',wt C write(0,*)'xa1',xa1 C write(0,*)'xb1',xb1 C write(0,*)'va1',pa1 C write(0,*)'vb1',pb1 C write(0,*)'xref',xref C C write(0,*)'VREF',vref C C write(0,*)'nt_inc=',nt_inc C write(0,*)'nt_ref=',nt_ref C write(0,*)'w_inc=',w_inc C write(0,*)'w_ref=',w_ref do it=1,nt_inc c { rtraj_match(1,it) = * rtraj_inc(1,it,1)*(1.-w_inc) + * rtraj_inc(1,it,2)*( w_inc) if(zmap .ne. 1) then rtraj_match(2,it) = * rtraj_inc(2,it,1)*(1.-w_inc) + * rtraj_inc(2,it,2)*( w_inc) else rtraj_match(2,it) = * rtraj_zmap_inc(2,it,1)*(1.-w_inc) + * rtraj_zmap_inc(2,it,2)*( w_inc) end if C rpar_match(1,it) = C * rpar_inc(1,it,1)*(1.-w_inc) + C * rpar_inc(1,it,2)*( w_inc) CC rpar_match(2,it) = C * rpar_inc(2,it,1)*(1.-w_inc) + C * rpar_inc(2,it,2)*( w_inc) rmod_match(1,it) = * rmod_inc(1,it,1)*(1.-w_inc) + * rmod_inc(1,it,2)*( w_inc) if(ray_in_tau) then rmod_match(2,it) = * rmod_inc(2,it,1)*(1.-w_inc) + * rmod_inc(2,it,2)*( w_inc) else rmod_match(2,it) = 0. end if C dt=sqrt(.25*((rtraj_match(2,it)-rtraj_match(2,it-1))- C * rmod_match(2,it)* C * (rtraj_match(1,it)-rtraj_match(1,it-1)))**2 + C * (1./rmod_match(1,it)* C * (rtraj_match(1,it)-rtraj_match(1,it-1)))**2 ) C write(0,*)'dt=',dt c } end do rtraj_match(1,nt_inc+1) = xref rtraj_match(2,nt_inc+1) = vref C! rpar_match(1,nt_inc+1) = pxref C! rpar_match(2,nt_inc+1) = pvref rpar_match(1,1) = pxref rpar_match(2,1) = pvref rmod_match(1,nt_inc+1) = velref if(ray_in_tau) then rmod_match(2,nt_inc+1) = sigref else rmod_match(2,nt_inc+1) = 0. end if rpar_match(1,2)=-rpar_ref(1,nt_ref,1)*(1.-w_ref) * -rpar_ref(1,nt_ref,2)*w_ref rpar_match(2,2)=-rpar_ref(2,nt_ref,1)*(1.-w_ref) * -rpar_ref(2,nt_ref,2)*w_ref do it=1,nt_ref c { rtraj_match(1,nt_inc+nt_ref+2-it) = * rtraj_ref(1,it,1)*(1.-w_ref) + * rtraj_ref(1,it,2)*( w_ref) if(zmap .ne. 1) then rtraj_match(2,nt_inc+nt_ref+2-it) = * rtraj_ref(2,it,1)*(1.-w_ref) + * rtraj_ref(2,it,2)*( w_ref) else rtraj_match(2,nt_inc+nt_ref+2-it) = * rtraj_zmap_ref(2,it,1)*(1.-w_ref) + * rtraj_zmap_ref(2,it,2)*( w_ref) end if C rpar_match(1,nt_inc+nt_ref+2-it) = C * -rpar_ref(1,it,1)*(1.-w_ref) - C * rpar_ref(1,it,2)*( w_ref) C rpar_match(2,nt_inc+nt_ref+2-it) = C * -rpar_ref(2,it,1)*(1.-w_ref) - C * rpar_ref(2,it,2)*( w_ref) rmod_match(1,nt_inc+nt_ref+2-it) = * rmod_ref(1,it,1)*(1.-w_ref) + * rmod_ref(1,it,2)*( w_ref) if(ray_in_tau) then rmod_match(2,nt_inc+nt_ref+2-it) = * rmod_ref(2,it,1)*(1.-w_ref) + * rmod_ref(2,it,2)*( w_ref) else rmod_match(2,nt_inc+nt_ref+2-it) = 0. end if c } end do do it=nt_inc+nt_ref+2,2*maxt c { rtraj_match(1,it) = rtraj_match (1,nt_inc+nt_ref+1) rtraj_match(2,it) = rtraj_match (2,nt_inc+nt_ref+1) C rpar_match(1,it) = rpar_match (1,nt_inc+nt_ref+1) C rpar_match(2,it) = rpar_match (2,nt_inc+nt_ref+1) rmod_match(1,it) = rmod_match (1,nt_inc+nt_ref+1) if(ray_in_tau) then rmod_match(2,it) = rmod_match (2,nt_inc+nt_ref+1) else rmod_match(2,it) = 0. end if c } end do if (srite('out',rtraj_match,2*nby_ray).ne.2*nby_ray) then call erexit("Problems writinf Output file") end if if (srite('Raypar_match',rpar_match,16) * .ne.16) then call erexit("Problems writinf Raypar_match file") end if if (srite('Raymod_match',rmod_match,2*nby_ray) * .ne.2*nby_ray) then call erexit("Problems writinf Raymod_match file") end if if (srite('Time_match',tmatch_tot,esize4) * .ne. esize4) then call erexit("Problems writing Time_match file") end if if (srite('Time_ref',tmatch_inc,esize4) * .ne. esize4) then call erexit("Problems writing Time_ref file") end if c } if (srite("location",location,12).ne.12) then call erexit("Trouble writing out location") end if if (srite("Pos_match",rtraj_match(1,nt_inc+1),8).ne.8) then call erexit("Trouble writing out Pos_match") end if end if c } end if c } end if c } end if c } end do c } end if c } end do if (.not. found) then c { bad=bad+1 write(0,*) "BAD RAY ALERT",bad,'sou_inc=',sou_inc+xsou,'sou_ref=',sou_ref+xsou if (srite('out',rtraj_match,2*nby_ray).ne.2*nby_ray) then call erexit("Problems writinf Output file") end if if (srite('Raypar_match',rpar_match,16).ne.16) then call erexit("Problems writinf Raypar_match file") end if if (srite('Raymod_match',rmod_match,2*nby_ray).ne.2*nby_ray) then call erexit("Problems writinf Raymod_match file") end if tmatch_tot=-1. tmatch_ref=-1. if (srite('Time_match',tmatch_tot,esize4) * .ne. esize4) then call erexit("Problems writing Time_match file") end if if (srite('Time_ref',tmatch_inc,esize4) * .ne. esize4) then call erexit("Problems writing Time_ref file") end if if (srite("location",location,12).ne.12) then call erexit("Trouble writing out location") end if if (srite("Pos_match",rtraj_match(1,nt_inc+1),8).ne.8) then call erexit("Trouble writing out Pos_match") end if end if C } end do C } end do end if c } c } end do if(want_multi .eq. 1) then write(0,*)'RECORD MULTIPLE ARRIVALS' call putch('n2','d',count_match+bad) call auxputch('n2','d',count_match+bad,'Raypar_match') call auxputch('n2','d',count_match+bad,'Raymod_match') call auxputch('n1','d',count_match+bad,'Time_match') call auxputch('n2','d',count_match+bad,'location') call auxputch('n1','d',count_match+bad,'Pos_match') else write(0,*)'FOUND MULTIPLE ARRIVALS BUT NOT RECORDES' end if end