#<
#Demultiple
#
#Usage:
#Demultiple.x  <in.H >out.H 
#
#Input Parameters:
#
#
#Output Parameters:
#
#
#
#Description:
#
#
#
#>
#%end of self-documentaiton
#-------------------------------------------------
#
#Author: Robert Clapp, ESMB 463, 7230253
#
#Date Created:Wed Apr 21 15:25:24 PDT 1999
#
#Purpose: 
#

program Demultiple{ 
	use sep3d_struct_mod
	use solver_mod+cgstep_mod+velsimp+nl_solver
  integer i1,i2,nv,niter,i3,niter_l,niter_nl
  real, allocatable :: vscan(:),cut(:),temp(:,:),tempv(:),tempd(:)
	real, pointer :: offset(:),data(:)
	real :: ov,dv,hcut
	integer :: nl,harold,mid,ia,ib,slp,auxin
	integer n2_c,n1_c,nsz,auxputch,putch,n3,n2
	real o1_c,d1_c,o2_c,d2_c,f1,f2,vel,fv,maxr
	integer pos1,pos2,iv,posv,sep3d_grab_key_index,ioff
	type(sep3d) :: in,out,mult
	character(len=1024) :: temp_ch

 	call init_3d()
  call init_sepf90()

	call init_sep3d("in",in,"INPUT")
	if(in%ndims>3) call erexit("expecting 3-D input (time, offset, cmp)")
	else if(in%ndims==3) n3=in%n(3)
	else{
		n3=1
		call sep3d_change_dims(in%tag,3,(/1,2,2/))
		call sep3d_grab_sep3d(in%tag,in)
	}


	if(0!=sep3d_grab_key_index(in%tag,"offset_x",ioff))
		call erexit("can not find offset key")

	call struct_init(in,out,"OUTPUT")
	call struct_init(in,mult,"OUTPUT")
	call sep3d_write_description("out",out)
	call sep3d_write_description("multiples",mult)
	call sep_get_header_format_tag("in",temp_ch)
	i1=putch("hff","s",temp_ch)
	i1=auxputch("hff","s",temp_ch,"multiples")
	call sep_get_grid_format_tag("in",temp_ch)
	i1=putch("gff","s",temp_ch)
	i1=auxputch("gff","s",temp_ch,"multiples")


	from par: integer nv,nl=0,harold=1,niter_l=5,mid=0,slp=5,niter_nl=10
	from par: real ov,dv,hcut

	if(-1==auxin("cut"))
		call erexit("must supply primary velocity function")

	from aux: cut integer n1:n1_c,n2:n2_c=1
	from aux: cut real o1:o1_c,d1:d1_c,o2:o2_c=0.,d2:d2_c=1.
	

	allocate(data(in%n(1)*in%n(2)),vscan(in%n(1)*nv),cut(n1_c*max(2,n2_c)),temp(in%n(1),in%n(2)))
	allocate(tempv(in%n(1)*nv),tempd(in%n(1)*in%n(2)),offset(in%n(2)))

	call sreed("cut",cut,n1_c*n2_c*4)
	if(n2_c==1) { n2_c=2; cut(n1_c+1:)=cut(1:n1_c)}


 	call ds_init(max_sample=10)


	call huber_cut(hcut)

	do i3=1,n3{
		call sep3d_grab_headers("in",in,nsz,fwind=(/0,i3-1/),nwind=(/in%n(2),1/))
		if(.not. sep3d_read_data("in",in,temp(:,:nsz)))
      call erexit("trouble reading data")
		if(size(offset)!=nsz){
			deallocate(data,tempd,offset)
			allocate(data(in%n(1)*nsz),tempd(in%n(1)*nsz),offset(nsz))
		}

		data=reshape(temp(:,1:nsz),(/in%n(1)*nsz/))
			

		call sep3d_grab_key_vals(in,ioff,offset)
		call velsimp_init(in%o(1),in%d(1),offset,ov,dv,in%n(1),nsz,nv)
	  to aux: vscan.H integer n1:in%n(1),n2:nv,n3:n3*2
		to aux: vscan.H real o1:in%o(1),o2:ov,d1:in%d(1),d2:dv
		vscan=0.

		if(i3==1) {
			call solver(velsimp_lop,cgstep,vscan,data,niter_l,verb=.true.)
			call nlcg_solver(velsimp_lop,huber,huber_init,vscan,data,niter_nl,
       x0=vscan,verb=.true.)
		}
		else{
				call solver(velsimp_lop,cgstep,vscan,data,niter_l,x0=vscan,verb=.true.)
		   call nlcg_solver(velsimp_lop,huber,huber_init,vscan,data,niter_nl,
         x0=vscan,verb=.true.)
		}

		call srite("vscan.H",vscan,size(vscan)*4)

		f2=(in%o(3)+in%d(3)*(i3-1)-o2_c)/d2_c+1.; pos2=f2; f2-=pos2
		if(pos2 < 1) {pos2=1;f2=0.}
		if(pos2 >= n2_c) {pos2=n2_c-1; f2=1.}
		tempv=0.
		do i1=1,in%n(1){
			f1=(in%o(1)+in%d(1)*(i1-1)-o1_c)/d1_c+1.; pos1=f1; f1-=pos1
			
			if(pos1 < 1) {pos1=1;f1=0.}
			if(pos1 >= n1_c) {pos1=n1_c-1; f1=1.}
		
			vel=(1.-f2)*(1.-f1)*cut(pos1+(pos2-1)*n1_c)
				+(f2)*(1.-f1)*cut(pos1+pos2*n1_c)
				+(f1)*(1.-f2)*cut(pos1+1+(pos2-1)*n1_c)
				+(f2)*f1*cut(pos1+1+pos2*n1_c)

		

			fv=(vel-ov)/dv+1.; posv=fv;fv-=posv
			if(posv < 1) { posv=1;fv=0.}
			if(posv >= nv) { posv=nv-1;fv=1.}

			
			do iv=1,posv
				tempv(i1+(iv-1)*in%n(1))=vscan(i1+(iv-1)*in%n(1))

			tempv(i1+posv*in%n(1))=fv*vscan(i1+posv*in%n(1))
		}
		call srite("vscan.H",tempv,size(vscan)*4)
		i1=velsimp_lop(.false.,.false.,tempv,tempd)
		call srite("multiples",tempd,size(tempd)*4)
		tempd=data-tempd
		call srite("out",tempd,size(tempd)*4)


		maxr=in%o(2)+in%d(2)*(in%n(2)-1)
		n2=2*maxr/in%d(2)+1
		to aux: full_both integer n1:in%n(1),n2:n2
		to aux: full_sig  integer n1:in%n(1),n2:n2
		to aux: full_both real o1:in%o(1),o2:-maxr
		to aux: full_sig  real o1:in%o(1),o2:-maxr
		to aux: full_both real d1:in%d(1),d2:in%d(2)
		to aux: full_sig  real d1:in%d(1),d2:in%d(2)
		deallocate(offset,tempd)
		allocate(offset(n2))
		do i2=1,n2
			offset(i2)=-maxr+(i2-1)*in%d(2)
		call velsimp_init(in%o(1),in%d(1),offset,ov,dv,in%n(1),n2,nv)
		allocate(tempd(in%n(1)*n2))
		i1=velsimp_lop(.false.,.false.,vscan,tempd)
		call srite("full_both",tempd,size(tempd)*4)
		i1=velsimp_lop(.false.,.false.,tempv,tempd)
		call srite("full_sig",tempd,size(tempd)*4)
		
	



	}
}