program PSpike_
  !! Paul Sava (paul@sep.stanford.edu)
  !! 09/21/99

  use sep
  use util_

  implicit none
  logical     :: verb

  type(myaxis):: a1,a2,a3,a4,a5,a6
  integer     :: n1,n2,n3,n4,n5,n6
  real        :: o1,o2,o3,o4,o5,o6
  real        :: d1,d2,d3,d4,d5,d6
  integer     :: i1,i2,i3,i4,i5,i6
  integer     :: nsp,esize,i,lk

  integer,dimension(:),allocatable :: k1,k2,k3,k4,k5,k6
  integer,dimension(:),allocatable :: l1,l2,l3,l4,l5,l6
  real,   dimension(:),allocatable :: mag, data

  call sep_init()

  call from_param("verb",verb,.false.)

  call from_param("n1",a1%n,1)
  call from_param("n2",a2%n,1)
  call from_param("n3",a3%n,1)
  call from_param("n4",a4%n,1)
  call from_param("n5",a5%n,1)
  call from_param("n6",a6%n,1)

  call from_param("o1",a1%o,0.)
  call from_param("o2",a2%o,0.)
  call from_param("o3",a3%o,0.)
  call from_param("o4",a4%o,0.)
  call from_param("o5",a5%o,0.)
  call from_param("o6",a6%o,0.)

  call from_param("d1",a1%d,1.)
  call from_param("d2",a2%d,1.)
  call from_param("d3",a3%d,1.)
  call from_param("d4",a4%d,1.)
  call from_param("d5",a5%d,1.)
  call from_param("d6",a6%d,1.)

  call houaxis(a1,1)
  call houaxis(a2,2)
  call houaxis(a3,3)
  call houaxis(a4,4)
  call houaxis(a5,5)
  call houaxis(a6,6)

  if(verb) then
     write(0,*) a1
     write(0,*) a2
     write(0,*) a3
     write(0,*) a4
     write(0,*) a5
     write(0,*) a6
  end if

  call from_param("esize",esize,4); call to_history("esize",esize)
  call from_param("nsp",nsp,0)

  if (nsp /=0 ) then
     allocate(k1(nsp),k2(nsp),k3(nsp),k4(nsp),k5(nsp),k6(nsp),mag(nsp))
     k1=1; call from_param("k1",k1,k1)!;if(maxval(k1)>a1%n) call seperr("k1>n1");
     k2=1; call from_param("k2",k2,k2)!;if(maxval(k2)>a2%n) call seperr("k2>n2");
     k3=1; call from_param("k3",k3,k3)!;if(maxval(k3)>a3%n) call seperr("k3>n3");
     k4=1; call from_param("k4",k4,k4)!;if(maxval(k4)>a4%n) call seperr("k4>n4");
     k5=1; call from_param("k5",k5,k5)!;if(maxval(k5)>a5%n) call seperr("k5>n5");
     k6=1; call from_param("k6",k6,k6)!;if(maxval(k6)>a6%n) call seperr("k6>n6");
     mag=1.; call from_param("mag",mag,mag)
     allocate(l1(nsp),l2(nsp),l3(nsp),l4(nsp),l5(nsp),l6(nsp))
     l1=a1%n; call from_param("l1",l1,l1)!;if(maxval(l1)>a1%n) call seperr("l1>n1");
     l2=a2%n; call from_param("l2",l2,l2)!;if(maxval(l2)>a2%n) call seperr("l2>n2");
     l3=a3%n; call from_param("l3",l3,l3)!;if(maxval(l3)>a3%n) call seperr("l3>n3");
     l4=a4%n; call from_param("l4",l4,l4)!;if(maxval(l4)>a4%n) call seperr("l4>n4");
     l5=a5%n; call from_param("l5",l5,l5)!;if(maxval(l5)>a5%n) call seperr("l5>n5");
     l6=a6%n; call from_param("l6",l6,l6)!;if(maxval(l6)>a6%n) call seperr("l6>n6");
     do i=1,nsp
        if(k1(i)>a1%n) k1(i)=a1%n
        if(l1(i)>a1%n) l1(i)=a1%n
        if(k1(i)< 1  ) k1(i)=1
        if(l1(i)< 1  ) l1(i)=1

        if(k2(i)>a2%n) k2(i)=a2%n
        if(l2(i)>a2%n) l2(i)=a2%n
        if(k2(i)< 1  ) k2(i)=1
        if(l2(i)< 1  ) l2(i)=1

        if(k3(i)>a3%n) k3(i)=a3%n
        if(l3(i)>a3%n) l3(i)=a3%n
        if(k3(i)< 1  ) k3(i)=1
        if(l3(i)< 1  ) l3(i)=1

        if(k4(i)>a4%n) k4(i)=a4%n
        if(l4(i)>a4%n) l4(i)=a4%n
        if(k4(i)< 1  ) k4(i)=1
        if(l4(i)< 1  ) l4(i)=1

        if(k5(i)>a5%n) k5(i)=a5%n
        if(l5(i)>a5%n) l5(i)=a5%n
        if(k5(i)< 1  ) k5(i)=1
        if(l5(i)< 1  ) l5(i)=1

        if(k6(i)>a6%n) k6(i)=a6%n
        if(l6(i)>a6%n) l6(i)=a6%n
        if(k6(i)< 1  ) k6(i)=1
        if(l6(i)< 1  ) l6(i)=1
     end do

     do i=1,nsp
        if(l1(i)<k1(i) .and. (l1(i)/=0 .and. k1(i)/=0)) then
           lk=k1(i); k1(i)=l1(i); l1(i)=lk
        end if
        if(l2(i)<k2(i) .and. (l2(i)/=0 .and. k2(i)/=0)) then
           lk=k2(i); k2(i)=l2(i); l2(i)=lk
        end if
        if(l3(i)<k3(i) .and. (l3(i)/=0 .and. k3(i)/=0)) then
           lk=k3(i); k3(i)=l3(i); l3(i)=lk
        end if
        if(l4(i)<k4(i) .and. (l4(i)/=0 .and. k4(i)/=0)) then
           lk=k4(i); k4(i)=l4(i); l4(i)=lk
        end if
        if(l5(i)<k5(i) .and. (l5(i)/=0 .and. k5(i)/=0)) then
           lk=k5(i); k5(i)=l5(i); l5(i)=lk
        end if
        if(l6(i)<k6(i) .and. (l6(i)/=0 .and. k6(i)/=0)) then
           lk=k6(i); k6(i)=l6(i); l6(i)=lk
        end if
     end do
     do i=1,nsp
        if (l1(i) .eq. 0 .and. k1(i) /=0) l1(i)=k1(i)
        if (l2(i) .eq. 0 .and. k2(i) /=0) l2(i)=k2(i)
        if (l3(i) .eq. 0 .and. k3(i) /=0) l3(i)=k3(i)
        if (l4(i) .eq. 0 .and. k4(i) /=0) l4(i)=k4(i)
        if (l5(i) .eq. 0 .and. k5(i) /=0) l5(i)=k5(i)
        if (l6(i) .eq. 0 .and. k6(i) /=0) l6(i)=k6(i)
     end do

  end if

  if(verb) then
     do i=1,nsp
        write(0,*) "---------------------------------"
        write(0,*) mag(i)
        write(0,*) k1(i),k2(i),k3(i),k4(i),k5(i),k6(i)
        write(0,*) l1(i),l2(i),l3(i),l4(i),l5(i),l6(i)
        write(0,*) "---------------------------------"
     end do
  end if

  call sep_close()
  allocate(data(a1%n))

  la6 :do i6=1,a6%n
     la5 :do i5=1,a5%n
        la4 :do i4=1,a4%n
           la3 :do i3=1,a3%n
              la2 :do i2=1,a2%n

                 data = 0.0
                 lsp:do i=1,nsp

                    if( &
                    k1(i)/=-1 .and. & 
                    k2(i)/=-1 .and. &
                    k3(i)/=-1 .and. &
                    k4(i)/=-1 .and. &
                    k5(i)/=-1 .and. &
                    k6(i)/=-1) then

                       if                (k1(i) /= 0) then
                          if             (k2(i) .eq. 0 .or. (i2>=k2(i).and.i2<=l2(i)) ) then
                             if          (k3(i) .eq. 0 .or. (i3>=k3(i).and.i3<=l3(i)) ) then
                                if       (k4(i) .eq. 0 .or. (i4>=k4(i).and.i4<=l4(i)) ) then
                                   if    (k5(i) .eq. 0 .or. (i5>=k5(i).and.i5<=l5(i)) ) then
                                      if (k6(i) .eq. 0 .or. (i6>=k6(i).and.i6<=l6(i)) ) then
                                         data(k1(i):l1(i)) = data(k1(i):l1(i)) + mag(i)
                                      end if !6
                                   end if !5
                                end if !4
                             end if !3
                          end if !2
                       else !1
                          if             (k2(i).eq.0 .or. (i2>=k2(i).and.i2<=l2(i)) ) then
                             if          (k3(i).eq.0 .or. (i3>=k3(i).and.i3<=l3(i)) ) then
                                if       (k4(i).eq.0 .or. (i4>=k4(i).and.i4<=l4(i)) ) then
                                   if    (k5(i).eq.0 .or. (i5>=k5(i).and.i5<=l5(i)) ) then
                                      if (k6(i).eq.0 .or. (i6>=k6(i).and.i6<=l6(i)) ) then
                                         data = data + mag(i)
                                      end if !6
                                   end if !5
                                end if !4
                             end if !3
                          end if !2
                       end if !1

                    end if !no spike

                 end do lsp
                 call sep_write (data)

              end do la2
           end do la3
        end do la4
     end do la5
  end do la6
  deallocate(data)
  if (nsp /=0 ) deallocate(k1,k2,k3,k4,k5,k6,mag,l1,l2,l3,l4,l5,l6)
  call exit (0)
end program PSpike_










