#<
#operator_calc
#
#Usage:
#operator_calc.x  <in.H >out.H 
#
#Input Parameters:
#
#
#Output Parameters:
#
#
#
#Description:
#
#
#
#>
#%end of self-documentaiton
#-------------------------------------------------
#
#Author: Bob Clapp
#
#Date Created:Mon Nov 16 10:38:48 1998
#
#Purpose: 
#

program operator_calc{ 
  integer nt,tau,nvel(2),nback,iloc(2),nv(2),i1,i,it,b(2),e(2),i2
	integer ia,ib
	real :: ovel(2),dvel(2),dt,time,ref(2),dray(2),tempr(2),tempr1(2),tempr2(2)
	real :: ov(2),dv(2),angle,d1(2),d2(2),len,ray_par(2),extra,h
  real,dimension(:,:,:), allocatable :: operator,vel2
  real,dimension(:,:), allocatable :: ray,vel,ray_mod,dz,slow0

	from history: integer n1:nt
	from history: real d1:dt,time
	from par: integer n1:nvel(1),n2:nvel(2),tau=0
	from par: real o1:ovel(1),o2:ovel(2)
	from par: real d1:dvel(1),d2:dvel(2)
	from aux: vel integer n1:nv(1),n2:nv(2)
	from aux: vel real o1:ov(1),o2:ov(2)
	from aux: vel real d1:dv(1),d2:dv(2)
	to history: integer n1:nvel(1),n2:nvel(2),n3:3,esize:4
	to history: real o1:ovel(1),o2:ovel(2)
	to history: real d1:dvel(1),d2:dvel(2)

	allocate(operator(nvel(1),nvel(2),3),ray(2,nt),vel(nv(1),nv(2)),ray_mod(2,nt),dz(nvel(1),nvel(2)),vel2(2,nv(1),nv(2)),slow0(nvel(1),nvel(2)))

	call sreed("in",ray,size(ray)*4)
	if(tau==0)
	call sreed("vel",vel,size(vel)*4)
	else{
		call sreed("vel",vel2,size(vel)*8)
		vel=vel2(1,:,:)
	}

	ref=ray(:,nt)
	operator=0.

	do i2=1,nvel(2){
 		iloc(2)=(ovel(2)+(i2-1)*dvel(2)-ov(2))/dv(2)+1.5
		do i1=1,nvel(1){
			ib=max(floor((ovel(1)+dvel(1)*i1-ov(1))/dv(1)+1.),1)
			ia=min(ceiling((ovel(1)+dvel(1)*i1-ov(1))/dv(1)+1.),nv(1))
			do i=ib,ia
				slow0(i1,i2)+=1./vel(i,iloc(2))
			slow0(i1,i2)=slow0(i1,i2)/(ia-ib+1)
		}
	}


	do i1=1,nt{
		iloc(1)=(ray(2,i1)-ov(1))/dv(1)+1.
		iloc(2)=(ray(1,i1)-ov(2))/dv(2)+1.
		ray_mod(1,i1)=vel(iloc(1),iloc(2))
	}
	d1=ray(:,nt)-ray(:,nt-1)
	angle=atan2(d1(1),d1(2))
	ray_par=(/cos(angle),sin(angle)/)
	ray_par=ray_par/ray_mod(1,nt)
	dz=0.



		

		do i1=2,nt-1{
			dray=ray(:,i1)-ray(:,i1-1)
			if(tau==0) len=dt*(ray_mod(1,i1)+ray_mod(1,i1-1))/2.
			else len=dray(1)**2/(dt*(ray_mod(1,i1)+ray_mod(1,i1-1))/2)


			tempr1=(/(ray(2,i1-1)-ovel(1))/dvel(1)+1,(ray(1,i1-1)-ovel(2))/dvel(2)+1.5/)
			b=tempr1
			tempr2=(/(ray(2,i1)-ovel(1))/dvel(1)+1,(ray(1,i1)-ovel(2))/dvel(2)+1.5/)
			e=tempr2

			if(b(1)!=e(1)){
				if(b(2)!=e(2)){
					call erexit("not expecting this")
				}
				else{
					tempr(1)=(ceiling(tempr1(1))-tempr1(1))/(tempr2(1)-tempr1(1))
					operator(e(1),e(2),1)+=len*(1.-tempr(1))
					operator(b(1),b(2),1)+=len*tempr(1)
				}
			}
			else{
				if(b(2)!=e(2)){
					tempr(2)=(ceiling(tempr1(2))-tempr1(2))/(tempr2(2)-tempr1(2))
					operator(b(1),b(2),1)+=len*(1.-tempr(2))
					operator(e(1),e(2),1)+=len*tempr(2)
				}
				else{
					operator(b(1),b(2),1)+=len
				}
			}
					
					



			#SECOND HALF
			tempr1=(/(ray(2,i1-1)-ovel(1))/dvel(1)+1,((ref(1)-ray(1,i1-1)+ref(1))-ovel(2))/dvel(2)+1.5/)
			b=tempr1
			tempr2=(/(ray(2,i1)-ovel(1))/dvel(1)+1,((ref(1)-ray(1,i1)+ref(1))-ovel(2))/dvel(2)+1.5/)
			e=tempr2

			if(b(1)!=e(1)){
				if(b(2)!=e(2)){
					call erexit("not expecting this")
				}
				else{
					tempr(1)=(ceiling(tempr1(1))-tempr1(1))/(tempr2(1)-tempr1(1))
					operator(e(1),e(2),1)+=len*(1.-tempr(1))
					operator(b(1),b(2),1)+=len*tempr(1)
				}
			}
			else{
				if(b(2)!=e(2)){
					tempr(2)=(ceiling(tempr1(2))-tempr1(2))/(tempr2(2)-tempr1(2))
					operator(b(1),b(2),1)+=len*(1.-tempr(2))
					operator(e(1),e(2),1)+=len*tempr(2)
				}
				else{
					operator(b(1),b(2),1)+=len
				}
			}
					
					
		}



			tempr=ray(:,nt)
			tempr=(/(tempr(2)-ovel(1))/dvel(1)+1.,(tempr(1)-ovel(2))/dvel(2)+1./)
			iloc=floor(tempr)
			dray=ray(:,nt)-ray(:,nt-1)
			if(tau==0) len=mod(time,dt)*ray_mod(1,nt)
			else len=dray(1)**2/(mod(time,dt)*ray_mod(1,nt))
			operator(iloc(1),iloc(2),1)+=len

			angle=ray_mod(1,nt)*ray_par(1)


			if(tau==0){
				it=iloc(1)
				extra=ray(2,nt)/dvel(1)+1.
				extra=(extra-floor(extra))
				do i1=2,it
					dz(i1,iloc(2))+=2*angle*slow0(i1-1,iloc(2))
				dz(it+1,iloc(2))+=2*angle*extra*slow0(it,iloc(2))

				do i2=1,nvel(2){
					operator(1,i2,2)=-dz(2,i2)/slow0(1,i2) * dvel(1)
					do i1=3,nvel(1)
						operator(i1-1,i2,2)=-dz(i1,i2)/slow0(i1-1,i2) * dvel(1)
				}

			}
		write(0,*) sum(operator(:,:,1)),sum(operator(:,:,2))

	call srite("out",operator,nvel(1)*nvel(2)*8)
	call srite("out",sum(operator,3),nvel(1)*nvel(2)*4)


}