anchor.diff:            call reflect(pos, min, dim, alpha)
             call findval(pos, value, dim)
             call minmax(value, min, max)
  -
  +          call centroid(pos, min, max, dim)
  +
  +          !! Expansion if the minimum value becomes the maximum
             if (minsave.EQ.max) then
  -              write(*,*) "expanding..."
  +              write(*,*) "expanding...", dist
  +             
                 maxsave = max
                 xsave = pos(1, max)
                 ysave = pos(2, max)
                 call expand(pos, max, dim, gamma)
                 call findval(pos, value, dim)
                 call minmax(value, min, max)
  -
  +              call centroid(pos, min, max, dim)
  +
  +              !! Setting things back, if expansion is worse than 
  +              !! standard reflection
                 if (maxsave.NE.max) then
  -                  pos(1, max) = xsave
  -                  pos(2, max) = ysave
  +                  pos(1, maxsave) = xsave
  +                  pos(2, maxsave) = ysave
  +                  call findval(pos, value, dim)
  +                  call minmax(value, min, max)
                 end if
   
  +          !! Contract from old position if minima value is still minima value
             else if (minsave.EQ.min) then
  -              write(*,*) "contracting..."
  +              !write(*,*) "contracting...", dist
                 xsave = pos(1, min)
                 ysave = pos(2, min)
   
                 pos(1, min) = xsave
                 pos(2, min) = ysave
                 call contract(pos, min, dim, beta)
  -
  +              call findval(pos, value, dim)
  +              call minmax(value, min, max)
  +              call centroid(pos, min, max, dim)
  +
  +              !! if minima is still minima, contract everything towards max
                 if (minsave.EQ.min) then
  +                  !write(*,*) "Contract all...", dist, min, max
  +                  !write(*,*) pos(1,1), pos(1,2), pos(1,3), pos(1,4)
  +                  !write(*,*) pos(2,:)
  +                  !write(*,*) pos
                     call contractall(pos, max, dim)
  +                  call findval(pos, value, dim)
  +                  call minmax(value, min, max)
  +                  call centroid(pos, min, max, dim)
  +
                 end if
   
             end if
             i = i + 1
  -      end do
  +          dist = sqrt((pos(1,min)-pos(1,max)) * (pos(1,min)-pos(1,max)) &
  +                      + (pos(2,min)-pos(2,max)) * (pos(2,min)-pos(2,max)))
  +
  +      end do
  +
  +      write(*,*) pos(1,1), pos(1,2), pos(1,3), pos(1,4)
  +      write(*,*) pos(2,1), pos(2,2), pos(2,3), pos(2,4)
  +      write(*,*) min, max, dim
   
   end subroutine
   
  
positive.code:       interface
            subroutine contract(pos, min, dim, beta)
                real*8, intent(inout)  :: pos(:,:)
                integer, intent(in)    :: min, dim
                real*8, intent(in)     :: beta
            end subroutine contract
        end interface
  
  
negative.code: 
        !!call init_random_seed()
  
        do i = 1,dim -1
            call random_number(pos(1,i))
            pos(1,i) = (pos(1,i) - 0.5)* 10
            call random_number(pos(2,i))
            pos(2,i) = (pos(2,i) - 0.5) * 10
        end do
        
  end subroutine
  
  !! finds centroid position 
  subroutine centroid(pos, min, max, dim)
        implicit none
        real*8, intent(inout)  :: pos(:,:)
        integer, intent(in)    :: min, max, dim
        integer                :: i
        real*8                 :: xsum, ysum
  
        pos(1, dim) = 0
        pos(2, dim) = 0
        xsum = 0
        ysum = 0
  
        do i = 1,dim - 1
            if (i.NE.min) then
               xsum = xsum + pos(1, i)
               ysum = ysum + pos(2, i)
            end if
        end do
  
        pos(1, dim) = xsum / real(dim - 2)
        pos(2, dim) = ysum / real(dim - 2)
        
  end subroutine
  
  