// This program draws the normal curve or "bell curve" used in statistics.
// You can pass in alternate number of points (default is 1000) to generate
// curves for different values.

use statistics.frink

plotNormal[mean, sigma, steps, g is graphics] :=
{
   low = 1/steps             // Use rational numbers so that the exactly
   high = 1-low              // right number of points is plotted.
   println["low is $low"]
   minSigma = inversePhi[low, 8]
   println["minsigma is $minSigma"]
   maxSigma = inversePhi[high, 8]
   println["maxsigma is $maxSigma"]
   
   vscale = 8 sigma^2           // Found experimentally to look good.
   ceilingH = normalDensity[mean + sigma * maxSigma, mean, sigma]
   scaledCeilingH = ceilingH * vscale
   r = scaledCeilingH
   g.stroke[r/2]

   println["ceiling H is $ceilingH"]
   println["scaled ceiling is $scaledCeilingH"]
   
   g.color[0.5,0.5,0.5]
   g.line[mean + (minSigma * sigma), 0, mean + (maxSigma * sigma), 0]

   // Draw SD lines
   highestSD = floor[maxSigma]
   for s = -highestSD to highestSD
   {
      x = mean + (sigma * s)
      y = -normalDensity[x, mean, sigma] * vscale
      g.line[x, 0, x, y]
   }
   width = maxSigma - minSigma

   // This polyline is the normal curve.
   c = new polyline
   for s=minSigma to maxSigma+0.001 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]

   wheel = r/2
   first = true
   points = 0
   for phi = high to low step ((low-high)/(steps-1))
   {
//      s = now[]
      x = inversePhi[phi,100,15]

      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, -1/2 r, r, r]
         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, r, r]

      points = points+1
//      e = now[]
//      println["point $points, time is " + format[e-s,"ms",3]]
   }
   println["$points points plotted."]
}

g = new graphics
points = 1000

// You can pass in a number of points as the sole argument.
if length[ARGS] > 0
   points = eval[ARGS@0]

plotNormal[100, 15, points, g]
g.show[]
g.write["normal$points.svg", 1024, undef]
g.write["normal$points.png", 2000, undef]
g.write["normal$points.html", 800, undef]
//g.print[]