--Copyright The Numerical Algorithms Group Limited 1991.

)r arrows
clipValue :SF :=  6
clipFun(x) == min(max(x,-clipValue),clipValue)
realSteps :INT := 25
imagSteps :INT := 25
C := Complex SF
S := Segment SF
drawComplexVectorField :(C -> C, S, S) -> VIEW3D
drawComplexVectorField(f,realRange,imagRange) ==
  delReal := (hi realRange - lo realRange)/realSteps
  delImag := (hi imagRange - lo imagRange)/imagSteps
  sp := create3Space()$ThreeSpace SF
  real := lo realRange
  for i in 1..realSteps + 1 repeat
    imag := lo imagRange
    for j in 1..imagSteps + 1 repeat
       z := f complex(real, imag)
       arg := argument z
       len := clipFun sqrt norm z
       p1 := point[real, imag, 0.0@SF, arg]
       scaleLen := delReal * len
       p2 := point[p1.1 + scaleLen * cos(arg),
                   p1.2 + scaleLen * sin(arg), 0.0@SF, arg]
       arrow := makeArrow(p1, p2)
       for a in arrow repeat curve(sp, a)
       imag := imag + delImag
    real := real + delReal
  makeViewport3D(sp, "Complex Vector Field")$VIEW3D

drawComplex :(C->C, S, S) -> VIEW3D
drawComplex(f, realRange, imagRange) ==
  deltaReal :SF := (hi realRange - lo realRange)/realSteps
  deltaImag :SF := (hi imagRange - lo imagRange)/imagSteps
  llp:List List Point SF := []
  real :SF := lo realRange
  for i in 1..realSteps + 1 repeat
    imag :SF := lo imagRange
    lp := []$(List Point SF)
    for j in 1..imagSteps + 1 repeat
      z :COMPLEX SF := f(complex(real, imag))
      pt :Point SF := point[real, imag, clipFun sqrt norm z, argument z]
      lp := cons(pt, lp)
      imag := imag + deltaImag
    real := real + deltaReal
    llp := cons(reverse! lp, llp)
  llp := reverse! llp
  makeViewport3D(mesh(llp), "Complex Function")$VIEW3D
