Computational changes in surffunc.f

leafres/testarea/surffunc.f

@@ -113,7 +113,7 @@
       double precision beta(5),root,der_root,fmax,yinter,der_yinter,
      &agivenx,der_agivenx,funcval_agivenx,xmin,xmax,curvatmax,xcurvatmax
       double precision a,b,c,x0,y0,term,term1,term2,term3,step,
-     &deratx,der2atx
+     &deratx,der2atx,yvector(5),xvector(5),zvector(5)
 
       a=beta(1)
       b=beta(2)
@@ -178,16 +178,39 @@
         der_root=-9999.0d0
       endif
       fmax=a+y0
-      call surffunc(1,yinter,1,0.0d0,ndim,beta,term,0)
-      call surffunc(1,term,1,0.0d0,ndim,beta,der_yinter,1)
-      call surffunc(1,term,1,agivenx,ndim,beta,der_agivenx,1)
-      call surffunc(1,funcval_agivenx,1,agivenx,ndim,beta,term,0)
+      xvector(1)=0.0d0
+      call surffunc(1,yvector(1),1,xvector(1),ndim,beta,zvector(1),0)
+      yinter=yvector(1)
+      call surffunc(1,yvector(1),1,xvector(1),ndim,beta,zvector(1),1)
+      der_yinter=zvector(1)
+      xvector(1)=agivenx
+      call surffunc(1,yvector(1),1,xvector(1),ndim,beta,zvector(1),1)
+      der_agivenx=zvector(1)
+      call surffunc(1,yvector(1),1,xvector(1),ndim,beta,zvector(1),0)
+      funcval_agivenx=yvector(1)
 
       curvatmax=-9999.0d0
       xcurvatmax=-9999.0d0
       step=(xmax-xmin)/1000.0d0
-      do term=xmin,xmax,step
-        call surffunc(1,term1,1,term,ndim,beta,deratx,1)
+      if(step.le.0.0d0)return
+!      do term=xmin,xmax,step
+!        call surffunc(1,term1,1,term,ndim,beta,deratx,1)
+!        term2=dexp(-(term-x0)/b)
+!        der2atx=-deratx/b+
+!     &(1.0d0+c)*deratx*deratx*((1.0d0+term2)**c)/(a*c)
+!        term3=dabs(der2atx/((1.0d0+deratx*deratx)**1.5d0))
+!        if(term3.gt.curvatmax)then
+!          curvatmax=term3
+!          xcurvatmax=term
+!        endif
+!      enddo
+
+      term=xmin
+      do while (term.le.xmax)
+        xvector(1)=term
+        call surffunc(1,yvector(1),1,xvector(1),ndim,beta,zvector(1),1)
+        term1=yvector(1)
+        deratx=zvector(1)
         term2=dexp(-(term-x0)/b)
         der2atx=-deratx/b+
      &(1.0d0+c)*deratx*deratx*((1.0d0+term2)**c)/(a*c)
@@ -196,7 +219,9 @@
           curvatmax=term3
           xcurvatmax=term
         endif
+        term=term+step
       enddo
+
       if(dabs(xcurvatmax-xmin).le.step.or.
      &dabs(xcurvatmax-xmax).le.step)then
         curvatmax=-9999.0d0