// This program draws the normal curve or "bell curve" used in statistics.
// It's a bit slow because calculating inverseErf for very high sigmas is
// quite slow.

use statistics.frink

plotNormal[mean, sigma, minSigma, maxSigma, g is graphics] :=
{
   g.color[0.5,0.5,0.5]
   vscale = 8 sigma^2           // Found experimentally to look good.
   g.line[mean + (minSigma * sigma), 0, mean + (maxSigma * sigma), 0]
   width = maxSigma - minSigma

   // This polyline is the normal curve.
   c = new polyline
   for s=minSigma to maxSigma step (width/100)
   {
      x = mean + (sigma * s)
      y = -normalDensity[x, mean, sigma] * vscale
      c.addPoint[x,y]
   }

   g.add[c]

   g.color[0,0,0]

   steps = 1000
   low = 1/1000                 // Use rational numbers so that the exactly
   high = 999/1000              // right number of points is plotted.
   wheel = 0
   first = true
   points = 0
   for phi = high to low step ((low-high)/(steps-1))
   {
      z = inversePhi[phi,8]
      x = mean + sigma * z

      n = normalDensity[x, mean, sigma]
      do
      {
         wheel = (wheel + 0.618034) mod 1
      } while wheel > n
      
      h = wheel
      if first
      {
         g.color[1,0,0]         // Draw the "you" circle in red.
         g.fillEllipseCenter[x, -.25, 1, 1]
         g.color[0,0,0]
         g.font["SansSerif", 4]
         g.text["You are here.", x, 7]

         g.line[x, 5, x, 1]     // Arrow body
         // Arrowhead
         p=new filledPolygon
         p.addPoint[x,.65]
         p.addPoint[x+0.3,2.5]
         p.addPoint[x-0.3,2.5]
         g.add[p]
         
         first = false
      } else
         g.fillEllipseCenter[x, -h*vscale, 1, 1]

      points = points+1
   }
   println["$points points plotted."]
}

g = new graphics
plotNormal[100, 15, -3.0902, 3.0902, g]
g.show[]
g.write["normal.svg", 800, 600]
g.write["normal.png", 800, 600]