dimension rea(40)
      dimension parvert(maxpar*(maxpar+1)),refl(maxpar),ssd(maxpar+1),
     +           expa(maxpar),cont(maxpar),var((maxvar+1)*maxpoints),
     +           par(maxpar),x(maxvar+1),center(maxpar)
c
      sdconv  = .false.
      parconv = .true.
      npar    = maxpar
      nvar    = maxvar
      npoints = 0
c
c     initialize the random number generator
c
      seed    = 72767
      xxx     = rand(seed)
c
c     testing to see if limits are exceeded.
c
      if ((npar.gt.maxpar).or.(nvar.gt.(maxvar+1)).or.(npoints.gt.
     +        maxpoints)) then
          write (*,1001) npar,maxpar,nvar,maxvar,npoints,maxpoints
          stop
      end if
c
c     reading the initial values for the parameters :
c
      read(*,1002) line
      call freeread(line,rea,nrea)
      do i=1,npar
          parvert(i)=rea(i)
      end do
c
c     reading all the values for the variables :
c
      ymean=0.
      ivar=0
   10 read(*,1002,end=20) line
          call freeread(line,rea,nrea)
          if (nrea.ge.(nvar+1)) then
             npoints=npoints+1
             do j=1,(nvar+1)
                ivar=ivar+1
                var(ivar)=rea(j)
                if (j.eq.(nvar+1)) ymean=ymean+var(nvar+1)
             end do
          endif
      goto 10
   20 continue
      ymean=ymean/npoints
c
c     generating the initial simplex :
c
c             first two vertex at +10% of first par and +/- 10% of all
c                            others.
c
      lastpar=npar
      do i1=1,npar
           i=i1*npar
           do j=1,npar
                parvert(i+j)=parvert(j)*(0.8+0.4*rand(0))
           enddo
      enddo
      lastpar=npar*(npar+1)
c
c     now that we have the initial vertex ready calculate the SSD for
c     each. first loop over all vertex :
c
      ncycles=-1
      ibest=1
      iworst=1
  100 ncycles=ncycles+1
         do i=1,(npar+1)
              ssd(i)=0.
              ifirstpar=(i-1)*npar
c
c             transfer all parameters from PARVERT to PAR :
c
              do j=1,npar
                    par(j)=parvert(ifirstpar+j)
              end do
c          
c             now call stddev for each vertex
c
              call stddev(par,var,ssd(i),x)
         end do
c
c        finding the best and worst vertex :
c
         do i=1,npar+1
             if (ssd(i).lt.ssd(ibest)) then
                  ibest=i
               else 
                  if (ssd(i).gt.ssd(iworst)) iworst=i
             end if
         end do
      ifirstbest=(ibest-1)*npar
      ifirstworst=(iworst-1)*npar
c
c     at this point IWORST contains the index of the worst vertex,
c     and IBEST that of the best.
c     test for convergence in any cycle other than the first.
c
      if (ncycles.ne.0) then
c
c        first condition is in SSD of best, RMS.LE.0.01% of YMEAN
c
         cond=abs(ymean/10000.)
         rms=sqrt(ssd(ibest)/npoints)
         if (rms.le.cond) sdconv=.true.
c
c        also stop if best and worst vertex are within 0.001% of each other
c
         parconv=.true.
         do i=1,npar
            i1=i*npar
            do j=1,npar
               if (parconv) then
                  test=parvert(i1+j)-parvert(j)
                  test=abs(100000.*test/parvert(j))
                  if (test.gt.1.0) parconv=.false.
               endif
            end do
         end do
c
         if ((parconv).or.(sdconv)) goto 999
      endif
c
c     reflect the worst vertex and obtain its SSD :
c
      ifirstworst=(iworst-1)*npar
      ifirstbest=(ibest-1)*npar
      do j=1,npar
           center(j)=0.
           do i=1,(npar+1)
                ifirstpar=(i-1)*npar
                if (i.ne.iworst) center(j)=center(j)+
     +                           parvert(ifirstpar+j)
           end do
           center(j)=parvert(ifirstworst+j)-(center(j)/(npar))
           refl(j)=parvert(ifirstworst+j)-(2.0*center(j))
      end do
      call stddev(refl,var,sdr,x)
c
c     compare the reflexed vertex with the best :
c
      if (sdr.le.ssd(ibest)) then
c
c            test expanded
c
           do i=1,npar
                expa(i)=refl(i)-center(i)
           end do
           call stddev(expa,var,sde,x)
           if (sde.le.sdr) then
c
c                here if expanded better than reflected better than best
c
                 ssd(iworst)=sde
                 do i=1,npar
                      parvert(ifirstworst+i)=expa(i)
                 end do
                 goto 100            
              else
c
c                here if expanded not better than reflected but reflected
c                better than best
c
                 ssd(iworst)=sdr
                 do i=1,npar
                       parvert(ifirstworst+i)=refl(i)
                 end do
                 goto 100
            end if
      endif
c
c     if reflected not better than best :
c
      if (sdr.lt.ssd(iworst)) then
c
c          if better than worst :
c
           ssd(iworst)=sdr
           do i=1,npar
              parvert(ifirstworst+i)=refl(i)
           end do
           goto 100
        else
c
c          if worse than worst try contracted
c        
           do i=1,npar
              cont(i)=parvert(ifirstworst+i)-(center(i)/2.)
           end do
           call stddev(cont,var,sdc,x)
           if (sdc.le.ssd(iworst)) then
c
c                 here if contracted better than worst
c
                  ssd(iworst)=sdc
                  do i=1,npar
                     parvert(ifirstworst+i)=cont(i)
                  end do
                  goto 100
               else
c
c                 no improvement, try shrinking
c
                  do i=1,(npar+1)
                  if (i.ne.ibest) then
                       ifp=(i-1)*npar
                       do j=1,npar
                          temp=parvert(ifp+j)
                          best=parvert(ifirstbest+j)
                          parvert(ifp+j)=temp-(0.5*(temp-best))
                       end do
                  end if
                  end do
                  goto 100
           end if
      end if
  999 write (*,1003) ncycles
      if (sdconv) write (*,1004)
      if (parconv) write (*,1005)
      write (*,'(1h )')
      do j=1,npar
          write (*,1006) j, parvert(ifirstbest+j)
      end do
      write (*,1007) rms,ssd(ibest)/npoints
c
c     to build the format for table :
c
      nvp3=nvar+3
      if (nvp3.gt.9) write (mf,'(i2)') nvp3
      if (nvp3.le.9) write (mf,'(i1)') nvp3
      wf='(1x,'//mf//'(d9.3,1x))'
c
c     to build the heading :
c
      do i=1,132
         hdng(i:i)=' '
      end do
      ihdng=2
      do i=1,nvar
          hdng(ihdng+3:ihdng+4)='x('
          write (mf,'(i2)') i
          hdng(ihdng+5:ihdng+6)=mf(1:2)
          hdng(ihdng+7:ihdng+7)=')'
          ihdng=ihdng+10
      end do
      hdng(ihdng+5:ihdng+5)='y'
      ihdng=ihdng+10
      hdng(ihdng+3:ihdng+7)='ycalc'
      ihdng=ihdng+10
      hdng(ihdng+4:ihdng+6)='dif'
      ihdng=ihdng+6
      write (*,'(a)') hdng(1:ihdng)
      write (*,'(1h )')
c
      reply='y'
      if ((reply(1:1).eq.'Y').or.(reply(1:1).eq.'y')) then
           do i=1,npar
                par(i)=parvert(ifirstbest+i)
           end do
           do i=1,npoints
              ifirstpoint=(i-1)*(nvar+1)
              do j=1,(nvar+1)
                 x(j)=var(ifirstpoint+j)
              end do
              yc=ycalc(par,x)
              dy=x(nvar+1)-yc
              write (*,fmt=wf) (x(j),j=1,(nvar+1)),yc,dy
           end do
      end if
      stop
 1001 format (/,5x,'npar=',i3,' maxpar=',i3,';  nvar=',i3,' maxvar=',
     +        i3,';   npoints=',i3,' maxpoints=',i3)
 1002 format(a256)
 1003 format (/,5x,'Convergence reached after',I4,' cycles',/)
 1004 format (10x,'RMS less than 0.1% of mean Y')
 1005 format (10x,'All the Vertex are within 0.01%')
 1006 format (5x,'a(',i2,') = ',D13.5)
 1007 format (/,10x,'RMS = ',D13.5,'  Std Dev = ',D13.5,///)
      end
c---------------------------------------------------------------------
      subroutine stddev(par,var,sd,x)
c---------------------------------------------------------------------
c  
c     It will loop over all points and calculate the summation of the
c     standard deviation for a given set of parameters.
c
      implicit double precision(a-h,o-z)
      common /nums/npar,nvar,npoints
      dimension par(2),var(2),x(2)
      sd=0.
c
c     looping over all points :
c
      do i=1,npoints
         ifirstv=(i-1)*(nvar+1)
c
c        transfering the variables from VAR to X :
c
         do j=1,(nvar+1)
              x(j)=var(ifirstv+j)
         end do
c
c        now calling the YCALC function to evaluate Y :
c
         yc=ycalc(par,x)
         dy=yc-x(nvar+1)
         sd=sd+(dy*dy)
      end do
      return
      end
c---------------------------------------------------------------------
      double precision function ycalc(a,x)
c---------------------------------------------------------------------
c
c     Evaluates YCALC for a given set of parameters and a set of
c     variables.
c
      implicit double precision(a-h,o-z)
      dimension a(*),x(*)
c
c     the definition of Y should be inserted here :
c