anchor.diff:          part1 = 1./(1.+ex(+.511*x/t9mev))
           part2 = 1./(1.+ex(-(x+2.531)*(.511/tnmev)-xi(1)))
   C       func2 = cnorm*(x+b*.511)*(x+2.531)**2*(x**2-1)**.5*part1*part2
  +      IF (x.gt.(2.531) THEN
  +        func2 = cnorm*(x-b*.511)*(x+2.531)**2*(x**2-1)**.5*part1*part2
  +      ELSE
           func2 = cnorm*(x+b*.511)*(x+2.531)**2*(x**2-1)**.5*part1*part2
  +      END IF
         END IF
         RETURN
         end
  
positive.code:         part1 = 1./(1.+ex(-.511*x/t9mev))
          part2 = 1./(1.+ex(+(x-2.531)*(.511/tnmev)-xi(1)))
  C       func1 = cnorm*x*(x-2.531)**2*(x**2-1)**.5*part1*part2
          func1 = cnorm*(x+b*.511)*(x-2.531)**2*(x**2-1)**.5*part1*part2
        END IF
        RETURN
        end
  
negative.code:         part1 = 1./(1.+ex(+.511*x/t9mev))
          part2 = 1./(1.+ex(-(x-2.531)*(.511/tnmev)+xi(1)))
  C       func4 = cnorm*x*(x-2.531)**2*(x**2-1)**.5*part1*part2
          func4 = cnorm*(X-b*.511)*(x-2.531)**2*(x**2-1)**.5*part1*part2
        END IF
        RETURN
        end
  