program da3d(da, da3dp, scene, output);
(* da3d: diana da file to 3d graphics

  Tom Schneider
  NCI/FCRDC Bldg 469. Room 144
  P.O. Box B
  Frederick, MD  21702-1201
  (301) 846-5581 (-5532 for messages)
  network address: toms@ncifcrf.gov

  National Cancer Institute
  Laboratory of Mathematical Biology

 *)

label 1; (* end of program *)

const
(* begin module version *)
version = 1.23; (* of da3d.p 1995 Nov 13
origin 1991 December 11 *)
(* end module version *)

(* begin module describe.da3d *)
(*
name
   da3d: diana da file to 3d graphics

synopsis
   da3d(da: in, da3dp: in, scene: out, output: out)

files
   da: output of the diana program; position to poistion correlations
   da3dp: parameter file to control scene.
           horizontal: shift of graph horizontally (in cm)

           vertical: shift of graph vertically (in cm)

           xlocation: location of viewer in bases

           ylocation: location of viewer in bases

           zlocation: location of viewer in bits

           magnify: magnification factor for whole scene, 1 = no change.

           xmagnify: magnification factor for x axis only.

           ymagnify: magnification factor for y axis only.

           zmagnify: magnification factor for z axis only.

           datacolumn: (integer) column of da which to use for the graph.

           dodiagonal: if the first character is 't', the diagonal
             values will be included.

   scene: 3D scene of the da data according to da3dp.  Result is in
      PostScript
   output: messages to the user

description
   Show the position to position correlation data in three dimensions.

examples

documentation

see also
   diana.p

author
   Thomas Dana Schneider

bugs

history
   [1995 Nov 13] Put lines to previous data trace so that it looks like a set
   of squares rather than just lines.

technical notes

*)
(* end module describe.da3d *)

(* begin module da3d.const *)
   mingrid = -100; (* minimum of side of the data grid *)
   maxgrid = +100; (* maxmum  of side of the data grid *)
(* end module da3d.const *)

(* begin module interact.const *)
      maxstring = 150; (* the maximum string *) 
(* end module interact.const version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.const *)
   pi = 3.14159265354; (* circumference divided by diameter of circle *)
   picfield = 8; (* width of numbers printed to the file *)
   picwidth = 5; (* number of decimal places for numbers *)
   charwidth = 0.0625; (* the width of characters in inches (ie, inches/char)
                        this allows centering of strings. *)
   (* note: for the Times-Roman  font, 0.0625 is a good value.
            for the Courier-Bold font, 0.08   is a good value. *)
   dotfactor = 0.00625; (* the size of dots *)
   defscale = 81; (* default scale factor.  coordinate units per inch *)
(* end module pic.const version = 2.63; (@ of dops.p 1991 November 2 *)

type

(* begin module pic.3d.type *)
(* these types are used by the three dimensional graphics routines *)

threevector = array[1..3] of real; (* a point in 3 space *)

tbtarray = array[1..3,1..3] of real; (* a three by three array *)

screen = record (* define a screen for viewing a 3d object *)
   a: threevector; (* center of screen *)
   b: threevector; (* screen x coordinate direction *)
   c: threevector; (* screen y coordinate direction *)
   v: threevector; (* the position of the viewer *)
   g: threevector; (* gaze: viewing direction *)
   smag: real; (* the magnification factor for the screen *)
   range: real; (* 1/smag; the half width of the screen *)
end;
(* end module pic.3d.type version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module interact.type *)
      string = record (* a string of characters *)  
         letters: array[1..maxstring] of char; (* the letters in the string *)
         length: integer; (* the number of characters in the string *)
         current: integer; (* the letter we are working on *) 
      end;  
(* end module interact.type version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module da3d.type *)
{   datagrid = array[mingrid,maxgrid..mingrid..maxgrid] of real;}
   datagrid = array[mingrid..maxgrid, mingrid..maxgrid] of real;
(* end module da3d.type *)

var
   da,   (* input file *)
   da3dp,  (* parameters *)
   scene  (* PostScript output *)
      : text; (* a file used by this program *)

(* begin module pic.var *)
   inpicture: boolean; (* true if we are drawing the picture,
                          ie, startpic has been called *)
   picxglobal, picyglobal: real; (* absolute location in the graph *)
   pictolerance: real; (* 10 raised to the picwidth,
      to detect values close to zero *)

   scale: real; (* scale factor.  graphic coordinate units per inch *)

(* NONSTANDARD for efficient use of postscript, keep track of
whether there is a current path *)
   inpath: boolean;

(* NONSTANDARD keep track of number of segments drawn so that
they can be stroked.  This (probably) solves the problem of the
Apple printer dying because it can't handle the data. *)
   segments: integer;

   xsideold, ysideold: real; (* current size of a rectangle.  see rectsize *)
(* end module pic.var version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module halt *)
procedure halt; 
(* stop the program.  the procedure performs a goto to the end of the   
   program.  you must have a label: 
      label 1;
   declared, and also the end of the program must have this label:  
      1: end. 
   examples are in the module libraries.
   this is the only goto in the delila system. *) 
begin 
      writeln(output,' program halt.'); 
      goto 1
end;  
(* end module halt version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module interact.clearstring *)
procedure clearstring(var ribbon: string);  
(* empty the string *)
var   index: integer; (* to the ribbon *) 
begin (* clearstring *) 
      with ribbon do begin  
         for index := 1 to maxstring do letters[index] := ' ';
         length := 0; 
         current := 0;
      end 
end; (* clearstring *)
(* end module interact.clearstring version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.functions *)
(* ********************************************************************** *)
(* begin module pic.await *)
procedure await;
(* Wait for user to type a carriage return.  the routine assumes that
there is a global file called input. *)
begin
(* the infinite way: 
    writeln(output);
    writeln(output,'*********************************');
    writeln(output,'* Use control-c to kill program *');
    writeln(output,'*********************************');
    while true do begin end;*)
end;
(* end module pic.await version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.startpic *)
procedure startpic(var afile:text; setscale,x,y: real; thefont: char);
(* open the graphics field, with the given scale, and at (x,y)
in that scale.  scale is in device coordinates per inch.
The font is chosen with thefont; t = Times-Roman, c = Courier-Bold *)
(* start pic output to file afile, set the globals *)
(* NONSTANDARD *)
(* this is the actual "world" coordinates used: *)
(*                        xmin,  xmax,  ymin, ymax *)
(*    ns; if (setwindow(-5.0/scale, +5.0/scale, -5.0/scale, +5.0/scale)*)
begin
   writeln(afile,'gsave'); (* save the current graphics state *)

   writeln(afile,'initgraphics'); (* make sure the printer
           is ready to print, without this, sometimes an Apple laserwriter
           will print the graph upside down, tiny and backwards! *)

   scale := setscale; (* set the global scale *)
   writeln(afile,'clear erasepage'); (* clean residue from before *)

   case thefont of
      'c': begin
           writeln(afile,'/Courier-Bold findfont'); (* locate the font *)
           writeln(afile,10:1,' scalefont'); (* set the font size in points*)
      end;
      't': begin
           writeln(afile,'/Times-Roman findfont'); (* locate the font *)
           writeln(afile,12:1,' scalefont'); (* set the font size in points*)
      end;
   end;

   writeln(afile,'setfont'); (* put the font into the current font *)

(* If the following statement is done then it will work on the sun, but will
kill the applewriter!!!! 

Sun's non-standard PostScript extension, setlinewidth has default 0,
as stated in the Read This First and the NeWS Manual.
This draws very quickly with 1 bit wide lines.
If '1 setlinequality' is not done, then one cannot set the width of lines.
So to use PostScript on the screen, I must first do
'1 setlinequality'.  However, if I send this code to the Applewriter,
it kills PostScript on the Applewriter and I get no output whatsoever!
(It took me several hours to figure this out, since once PostScript is killed
on the Applewriter, the NEXT output is also smashed and I had to figure that
out also...)

So a standard PostScript program will not work correctly with the
default.  "Correcting" the PostScript program so that it works on the
Sun means that it BOMBS on the Applewriter.

The default for setlinequality should be '1 setlinequality' so that
the same PostScript code can be used both on the Sun and with other devices.
If you want speed, use the nonstandard form.
An alternative is to redefine setlinequality so that '0 setlinequality'
does give correct results with standard PostScript.

Please review this, Randy.  I think that Sun should fix it.
   writeln(afile,'1 setlinequality');  makes lines at least 1 bit wide *)

   (* set the scale to inches
   writeln(afile,
                 scale:picfield:picwidth,' ',
                 scale:picfield:picwidth,' scale'); *)

   (* define some things in postscript *)
   (* doline allows less stuff to be put in the output file.
      it takes two numbers off the stack, copies them, draws a line
      to them as coordinates. *)
(* replaced by 'currentpoint translate'
   writeln(afile,'/doline { 2 copy lineto } def');
*)
   (* define a function that makes inches out of a number *)
(* do this all internally here, it's faster
   writeln(afile,'/i { ',scale:picfield:picwidth,' mul} def');
*)

   (* move to the start point on the page *)
   writeln(afile,
                     (x*scale):picfield:picwidth,
                 ' ',(y*scale):picfield:picwidth,
                 ' translate');

    writeln(afile);
    writeln(afile,'% Define functions so the text produced is smaller');
    writeln(afile,'/a {stroke newpath 0 0} def     % special for arc');
    writeln(afile,'/c {stroke 0 0 moveto} def      % current point');
    writeln(afile,'/f {findfont 10 scalefont setfont} def');
    writeln(afile,'   % to set fonts simply use the f function.  Example:');
    writeln(afile,'   %/Symbol f (\142) /Courier-Bold f (-galactosidase');
    writeln(afile,'/l {lineto} def');
    writeln(afile,'/m {moveto} def');
    writeln(afile,'/n {stroke newpath 0 0 moveto} def'); (* new segment *)
    writeln(afile,'/rl {rlineto} def');
    writeln(afile,'/rm {rmoveto} def');
    writeln(afile,'/s {newpath 0 0 moveto} def     % Start path '); 
    writeln(afile,'/t {currentpoint translate} def % translate ');
    writeln(afile,'/x {show} def                   % show teXt ');
    writeln(afile);

   (* start out the pathway *)
   inpath := false;
   (* start the number of segments written: *)
   segments := 0;

   (* now for the normal pic stuff: *)
   inpicture := true;
   picxglobal := 0.0;
   picyglobal := 0.0;
   pictolerance := trunc(exp(picwidth*ln(10))+0.5)
(*;writeln(output,'pictolerance = ',pictolerance:picfield:picwidth);*)
end;
(* end module pic.startpic version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.stoppic *)
procedure stoppic(var afile:text);
(* stop pic output to file afile *)
(* NONSTANDARD *)
begin
    if inpath then begin
       writeln(afile,'stroke');
       inpath := false
    end;

    writeln(afile,'showpage');
    writeln(afile,'grestore'); (* restore the current graphics state
                                  to what it was before the startpic *)

    await;

    inpicture := false;
end;
(* end module pic.stoppic version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.drawr *)
procedure drawr(var afile: text; dx,dy: real; visibility: char;
               spacing: real);
(* make a line to file afile by relative draw of dx,dy with visibility
  i invisible
  - dashed
  . dotted
  l line
with the dashes or dots separated by the spacing given
(this has no effect with invisible and line). *)
(* NONSTANDARD *)
var
   ddx,ddy: real; (* changes in dx and dy for dots and dashes *)
   dr: real; (* the hypotenuse, the distance actually drawn *)
   on: boolean; (* draw linesegment if true *)
   y: real; (* the variable for tracking dots and dashes *)
   r: integer; (* number of times to cycle for dots and dashes *)
   ss: real; (* precalculated value to make things a bit faster *)
   theta: real; (* angle of the line *)
procedure checkseg(var afile: text);
(* NONSTANDARD checks how many segments have been written, if
more than 'buffer', stroke them to the postscript page *)
const buffer = 10;
begin
   if segments >= buffer
   then begin
   (* New segment: writeln(afile,'stroke newpath 0 0 moveto'); *)
      writeln(afile,'n');
      segments := 0
   end
   else segments := segments + 1;
end;
begin (* drawr *)
   if not inpath then begin
      (* starts from current coordinates *)
      (* Start path: writeln(afile,'newpath 0 0 moveto'); *)
      writeln(afile,'s');
      inpath := true
   end
   else checkseg(afile);
(* checks
if not (visibility in ['l','i','.','-'])
then writeln(afile,'%YELLLLLL!!!',visibility,'!');
writeln(afile,'% ',visibility,' line');*)

   (* put these on the stack, they will always be used *)
   write(afile,    (dx*scale):picwidth:picfield,
               ' ',(dy*scale):picwidth:picfield);
   case visibility of
  'l','i': begin
              case visibility of
                 'i': write(afile,' m');
                 'l': write(afile,' l');
              end
           end;
  '.','-': begin (* make up our own dots and dashes *)
              writeln(afile); (* move away from the (dx,dy) on the stack *)
              if spacing <= 0.0 then begin
                 writeln(output,'drawr: spacing zero with . or - line');
                 halt
              end;
              if dx = 0.0
              then begin
                 ddx := 0.0; (* avoid division by zero *)
                 ddy := scale*spacing;
                 if dy < 0 then ddy := - ddy; (* this makes sure that
                    we draw lines straight down if that was the request *)
              end
              else begin
                 (* find out the angle of the slope, intentionally
                    lose the sign *)
                 theta := arctan(abs(dy/dx));

                 ddx := scale*spacing*cos(theta);
                 ddy := scale*spacing*sin(theta);

                 (* return the sign to the little buggers *)
                 if dx < 0 then ddx := -ddx;
                 if dy < 0 then ddy := -ddy;
              end;
              y := 0;
              case visibility of
                 '.': ss := scale*dotfactor;
                 '-': on := true;
              end;

              dr := sqrt(dx*dx+dy*dy);
              for r := 1 to round(dr/spacing) do begin
                 case visibility of
                 '-': begin
                         write(afile,    (ddx):picwidth:picfield,
                                     ' ',(ddy):picwidth:picfield);
                         if on
                         then writeln(afile,' rl')
                         else writeln(afile,' rm');
                         on := not on
                      end;
                 '.': begin
                         (* put out a dot like in dotr *)
                         write(afile,
                                    +ss:picwidth:picfield,' 0 rl');
                         write(afile,' ',
                                    -ss:picwidth:picfield,' 0 rl');
                         write(afile,' ',(ddx):picwidth:picfield,
                                     ' ',(ddy):picwidth:picfield);
                         writeln(afile,' rm');
                      end;
                 end
              end;
              (* let's make really sure we got there!! *)
              writeln(afile,' m'); (* pulled from the stack *)
           end;
   end;

   (* an elegant way to make postscript keep a global record is
   to translate the coordinates! *)
   (* writeln(afile,' currentpoint translate'); *)
   writeln(afile,' t');

   picxglobal := picxglobal + dx;
   picyglobal := picyglobal + dy;
end;
(* end module pic.drawr version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.mover *)
procedure mover(var afile: text; dx,dy: real);
(* move relative the amount (dx, dy). *)
begin
   drawr(afile,dx,dy,'i',0.0);
end;
(* end module pic.mover version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.liner *)
procedure liner(var afile: text; dx,dy: real);
(* draw a line the relative amount (dx, dy). *)
begin
   drawr(afile,dx,dy,'l',0.0);
end;
(* end module pic.liner version = 2.63; (@ of dops.p 1991 November 2 *)


(* begin module pic.drawa *)
procedure drawa(var afile: text; x,y: real; visibility: char;
               spacing: real);
(* make a line to file afile to absolute coordinate x,y with visibility
  i invisible
  - dashed
  . dotted
  l line
with the dashes or dots separated by the spacing given
(this has no effect with invisible and line). *)
var
   dx, dy: real; (* differences between current and desired
                    locations *)
begin
   dx := x - picxglobal;
   dy := y - picyglobal;

   drawr(afile,dx,dy,visibility,spacing)
end;
(* end module pic.drawa version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.movea *)
procedure movea(var afile: text; x,y: real);
(* move to absolute x and y *)
begin
   drawa(afile,x,y,'i',0.0);
end;
(* end module pic.movea version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.linea *)
procedure linea(var afile: text; x,y: real);
(* draw a line from current position to absolute x and y *)
begin
   drawa(afile,x,y,'l',0.0);
end;
(* end module pic.linea version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.graphstring *)
procedure graphstring(var tofile: text; var s: string; centered: boolean); 
(* graph the string s.  If it is recognized as a quoted string (surrounded
by double quotes), graph it without the quotes and center it.
Always center if centered is true.
Otherwise simply graph it.  if not in picture, just write it to output *)
(* NONSTANDARD *)
var i: integer; (* index to s, and temporary storage *)
    mv: real; (* holds amount to move, in plotting coordinates *)
    quoted: boolean; (* true if the string is quoted *)
begin
   if (inpicture and (s.length > 0))
   then with s do begin
      if length > 2
      then if (letters[1]='"') and (letters[length]='"')
           then quoted := true
           else quoted := false
      else quoted := false;

      (* override so quoted strings are always centered *)
      if quoted then centered := true;

      if centered then begin
         (* generate the calls to center the string.  Note: this must not
            be done for the pic program, which already centers *)
         if quoted then i := length-2
                   else i := length;
         mv := i*charwidth/2.0;
         mover(tofile,-mv,0.0);
      end;

      (* do the non-standard postscript: *)
      (* do postscript to complete pervious path *)
      (* set current point: writeln(tofile,'stroke 0 0 moveto'); *)
      writeln(tofile,'c');

      write(tofile,'('); (* begin postscript literal *)
      if quoted
      then for i := 2 to length-1 do write(tofile,letters[i])
      else for i := 1 to length   do write(tofile,letters[i]);
      write(tofile,')'); (* end postscript literal *)
      writeln(tofile,' x'); (* show the literal *)
      inpath := false; (* force new path from here *)

      if centered then begin
         (* restore to previous location *)
         mover(tofile,mv,0.0);
      end;
   end
   else begin
(* no output if not in picture 
      writestring(tofile,s);
      writeln(tofile) *)
   end
end;
(* end module pic.graphstring version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.stringinteger *)
procedure stringinteger(number: integer; var name: string;
                        width: integer; leadingzeros: boolean);
(* make the string from the number, start putting characters in
after the current length point. use width characters.
if leadingzeros is true, trail zeros before the number. *)
var
   bigdigit: integer; (* the location of the biggest digit *)
   dig: integer; (* number of digits in the number *)
   place: integer; (* place to write the next digit of the number *)
   sign: integer; (* the sign of the number *)
begin
   with name do begin
      if number < 0
      then begin
         sign := -1;
         length := length + 1; (* provide room for the sign!! *)
         number := -number;
         if leadingzeros then begin
            writeln(output,'WARNING: stringinteger: the sign of a negative',
                           ' number with leading zeros is lost');
         end
      end
      else sign := +1;

      (* log 10 of the number plus 1 is the number of digits in the number.
      On this sun computer ln(1000)/ln(10) is 2.9999, which when
      truncated gives 2, rather than the desired 3.  To avoid this
      kind of problem, 0.1 is added. *)
      if number > 9
      then dig := trunc(ln(number+0.1)/ln(10))+1
      else dig := 1;

      if dig > width then begin
         writeln(output,'stringinteger: number width too small');
         writeln(output,dig:1,' digit number (',number:1,')');
         writeln(output,'does not fit in ',width:1,' characters');
         halt
      end;
      if leadingzeros
      then bigdigit := length + 1 (* no sign if leading zeros *)
      else begin
         bigdigit := length + width - dig + 1;

         if (bigdigit <= length) and (sign < 0) then begin
            writeln(output,'stringinteger: no room for sign');
            halt
         end;
      end;
      if sign < 0 then letters[bigdigit-1] := '-';

      for place := length + width downto bigdigit do begin
         case (number mod 10) of
            0: letters[place] := '0';
            1: letters[place] := '1';
            2: letters[place] := '2';
            3: letters[place] := '3';
            4: letters[place] := '4';
            5: letters[place] := '5';
            6: letters[place] := '6';
            7: letters[place] := '7';
            8: letters[place] := '8';
            9: letters[place] := '9';
         end;
         number := number div 10;
      end;
      length := length + width;
   end
end;
(* end module pic.stringinteger version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.stringreal *)
procedure stringreal(number: real; var name: string;
                     width, decimal: integer);
(* make the string from the real number, start putting characters in
at the start point. use width characters and decimal characters
after the decimal place *)
   (* note that the rounding operation to get the digits below zero
      must be done first.  then the digits above zero can be lopped off.
      this makes 99.99 come out correctly to 100.0 (to 1 decimal place)
      otherwise, 99.99 -> 0.99 -> 1.0 (rounded) -> 10 (print with 1 decimal
      place), and stringinteger won't be happy about that. *)
var
   abovezero: integer; (* the number shifted above the decimal place, to 
      'decimal' positions (and rounded) *)
   shift: integer; (* power of ten used to shift a number around
      relative to the decimal point *)
   sign: integer; (* the sign of the number *)
   thedecimal: integer; (* integer version of the decimal part of the number *)
   theupper: integer; (* integer version of the upper part of the number *)
begin
   if number < 0 then sign := -1
                 else sign := +1;

   number := abs(number); (* make positive *)

   (* the amount to shift the number above zero *)
   shift := round(exp(decimal*ln(10))); (* amount to move above zero *)
   abovezero := round(number*shift); (* move above zero, round off *)
   theupper := trunc(abovezero/shift);
   thedecimal := abovezero - shift*theupper;

   (* create the actual real number *)
   (* before decimal point *)
   stringinteger(sign*theupper,name,width-decimal-1,false);
   with name do begin (* put in the decimal point *)
      length := length + 1;
      letters[length] := '.';
   end;
   stringinteger(thedecimal,name,decimal,true); (* after decimal point *)
end;
(* end module pic.stringreal version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.picnumber *)
procedure picnumber(var afile: text;
                    dx, dy, number: real; width, decimal: integer;
                    centered: boolean);
(* Supply graphic commands for a 'number' whose center is at the relative point
(dx, dy) from the current point, 'width' characters wide and 'decimal'
characters beyond the decimal point.
If the width is zero, no number is produced.
procedure stringnumber(number: integer; start: integer; var name: string);
the location after the call is the same as before the call.
The string is optionally centered *)
var
   name: string; (* the string to pack the number into for shipping out *)
begin
   if width > 0 then begin
      mover(afile,dx,dy);

      clearstring(name);

      if decimal>0
      then stringreal(number,name,width,decimal)
      else stringinteger(round(number),name,width,false);

      graphstring(afile, name, centered);
      mover(afile,-dx,-dy);
   end
end;
(* end module pic.picnumber version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.xtic *)
procedure xtic(var afile: text; length, dx, dy, number: real;
               width, decimal: integer);
(* produce a tic mark for the x axis of "length" long.
Supply a number whose center is at the relative point (dx, dy)
from the end to the tick, 'width' characters wide and 'decimal'
characters beyond the decimal point.
If the width is zero, no number is produced.
the location after the call is the same as before the call. *)
begin
   liner(afile,0.0,-length);
   picnumber(afile,dx,dy,number,width,decimal,true);
   mover(afile,0.0,length);
end;
(* end module pic.xtic version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.ytic *)
procedure ytic(var afile: text; length, dx, dy, number: real;
               width, decimal: integer);
(* produce a tic mark for the y axis of "length" long.
Supply a number whose center is at the relative point (dx, dy)
from the end to the tick, 'width' characters wide and 'decimal'
characters beyond the decimal point.
If the width is zero, no number is produced.
the location after the call is the same as before the call. *)
begin
   liner(afile,-length,0.0);
   picnumber(afile,dx,dy,number,width,decimal,true);
   mover(afile,length,0.0);
end;
(* end module pic.ytic version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.xaxis *)
procedure xaxis(var afile: text;
         axlength,fromtic,interval,totic: real;
         length, dx, dy: real;
                width, decimal: integer);
(* draw an x axis starting from the current position.
the length of the xaxis is axlength.
the axis is labeled with numbers starting with fromtic
at intervals given up to totic.
the remaining variables describe the form of the tic marks as in xtic.
If the width is zero, no number is produced.
the location after the call is the same as before the call. *)
var
   jump: real; (* the space to move on the graph between tic marks *)
   jumpdistance: real; (* the total jumps made.  this may not be
      a simple function of the input variables since they may
      not work out to an exact number of jumps *)
   tic: real; (* the numerical value of the tic label *)
   totici2: real; (* an extra half interval guarantees the last tic is placed *)
begin
   liner(afile,axlength,0.0);
   mover(afile,-axlength,0.0);
   if totic = fromtic then begin
      writeln(output,'xaxis: fromtic and totic cannot be equal');
      halt;
   end;
   if (axlength = 0.0) or (interval = 0.0) then begin
      writeln(output,'xaxis: neither axlength nor interval can be zero');
      halt;
   end;
   jump := axlength * interval / (totic - fromtic);
   jumpdistance := 0;

   (* note!  the extra amount of interval/2 makes sure that the last tic
      is placed, if the computer does not add quite right and overshoots
      by a tich on the value of tic.  This much makes it surely work! *)
   tic := fromtic;
   if interval > 0.0 then begin
      totici2 := totic + interval/2.0;
      while tic <= totici2 do begin
         xtic(afile,length,dx,dy,tic,width,decimal);
         tic := tic + interval;
         if tic <= totici2 then begin
            mover(afile,jump,0.0);
            jumpdistance := jumpdistance + jump;
         end
      end
   end
   else if interval < 0.0 then begin
      totici2 := totic - interval/2.0;
      while tic >= totici2 do begin
         xtic(afile,length,dx,dy,tic,width,decimal);
         tic := tic + interval;
         if tic >= totici2 then begin
            mover(afile,jump,0.0);
            jumpdistance := jumpdistance + jump
         end
      end
   end;

   mover(afile,-jumpdistance,0.0)
end;
(* end module pic.xaxis version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.yaxis *)
procedure yaxis(var afile: text;
         aylength,fromtic,interval,totic: real;
         length, dx, dy: real;
                width, decimal: integer);
(* draw a y axis starting from the current position.
the length of the yaxis is aylength.
the axis is labeled with numbers starting with fromtic
at intervals given up to totic.
the remaining variables describe the form of the tic marks as in ytic.
If the width is zero, no number is produced.
the location after the call is the same as before the call. *)
var
   half: real; (* half the jump interval.  By adding this to the while loops,
                  we assure that the very last tic gets done, and isn't lost
                  due to roundoff *)
   jump: real; (* the space to move on the graph between tic marks *)
   jumpdistance: real; (* the total jumps made.  this may not be
      a simple function of the input variables since they may
      not work out to an exact number of jumps *)
   tic: real; (* the numerical value of the tic label *)
begin
   liner(afile,0.0,aylength);
   mover(afile,0.0,-aylength);
   if totic = fromtic then begin
      writeln(output,'yaxis: fromtic and totic cannot be equal');
      halt;
   end;
   if (aylength = 0.0) or (interval = 0.0) then begin
      writeln(output,'yaxis: neither aylength nor interval can be zero');
      halt;
   end;
   jump := aylength * interval / (totic - fromtic);
   jumpdistance := 0;
   half := interval / 2.0;

   tic := fromtic;
   if interval > 0.0 then while tic <= totic+half do begin
      ytic(afile,length,dx,dy,tic,width,decimal);
      tic := tic + interval;
      if tic <= totic+half then begin
         mover(afile,0.0,jump);
         jumpdistance := jumpdistance + jump
      end
   end
   else if interval < 0.0 then while tic >= totic-half do begin
      ytic(afile,length,dx,dy,tic,width,decimal);
      tic := tic + interval;
      if tic >= totic-half then begin
         mover(afile,0.0,jump);
         jumpdistance := jumpdistance + jump
      end
   end;

   mover(afile,0.0,-jumpdistance)
end;
(* end module pic.yaxis version = 2.63; (@ of dops.p 1991 November 2 *)
(* ********************************************************************** *)
(* end module pic.functions version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.polrec *)
procedure polrec(r,theta: real; var x,y: real);
(* convert polar to rectangular coordinates,
theta is in radians *)
begin
   x := r*cos(theta);
   y := r*sin(theta)
end;
(* end module pic.polrec version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.3d.package *)
(* ********************************************************************** *)
(* begin module pic.3d.determinant *)
function determinant(a: tbtarray): real;
(* compute the determinant of a *)
begin
   determinant := +a[1,1] * (a[2,2]*a[3,3] - a[3,2]*a[2,3])
                  -a[1,2] * (a[2,1]*a[3,3] - a[3,1]*a[2,3])
                  +a[1,3] * (a[2,1]*a[3,2] - a[3,1]*a[2,2])
end;
(* end module pic.3d.determinant version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.3d.d32 *)
procedure d32(o, a, b, c, v: threevector;  var xloc,yloc: real);
(* convert from 3d to 2d.  the players are:
o: the coordinate of the object point to be converted to 2d
a,b,c: define the position of the window (screen):
a: center of screen
b: screen x coordinate direction
c: screen y coordinate direction
v: the position of the viewer
xloc,yloc: the resulting image vector in screen coordinates.
   The method of graphics is to project the object (o) toward the viewer
(v) and to determine the interception of this line with the screen
as defined by a,b and c.  the result is expressed in the coordinate system
of the screen, and so can be plotted on a 2d plotting device.
   When one works through the vector math, it turns out that to find
the screen coordinates requires solving a set of linear equations.
This is done using Cramer's rule and determinants. *)
var
   ov,oa: real; (* for partial calculation *)
   j: integer; (* index to the arrays *)
   d,x,y: tbtarray;
begin
   (* define the coefficients of the equations in d,x and y *)
   for j:=1 to 3 do begin
      ov := o[j]-v[j];
      d[j,1]:=b[j];
      d[j,2]:=c[j];
      d[j,3]:=ov;

      oa:=o[j]-a[j];
      x[j,1]:=oa;
      x[j,2]:=c[j];
      x[j,3]:=ov;

      y[j,1]:=b[j];
      y[j,2]:=oa;
      y[j,3]:=ov;
   end;

   (* use cramer's rule to find the solution *)
   xloc:=determinant(x)/determinant(d);
   yloc:=determinant(y)/determinant(d);
end;
(* end module pic.3d.d32 version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.3d.view *)
procedure view(v: threevector; var gaze: threevector; smag: real;
               var a,b,c: threevector);
(* this routine converts a viewing position (v) and a viewing
direction (gaze), into the a,b,c values of a vertically oriented screen
(ie, the screen is right side up).  a is the center of the screen,
b is the x axis, c is the y axis on the screen.  This saves the user
the trouble to make sure that b, c and the direction of viewing are
orthogonal.

one may magnify the view by making smag greater than one, or one may
shrink the view by making smag less than one.

if the viewing direction vector is not large enough,
then the program halts.

note: gaze is automatically converted to a unit vector. *)
var
   db: real; (* magnitude of db *)
   dgaze: real; (* magnitude of gaze *)
   j: integer; (* index to the arrays *)
begin
   (* first check out the gaze direction *)
   dgaze := sqrt(gaze[1]*gaze[1] + gaze[2]*gaze[2] + gaze[3]*gaze[3]);
   if smag = 0.0 then begin
      writeln(output,'screen magnitude cannot be zero');
      halt
   end;
   if dgaze <= 0.001 then begin
      writeln(output,'gaze magnitude (',dgaze:5:3,') is too small');
      halt
   end;

   (* make gaze a unit vector and set up the a vector as the
      viewing point plus the gaze vector *)
   for j := 1 to 3 do begin
      gaze[j] := gaze[j]/dgaze;
      a[j] := v[j] + gaze[j]
   end;

   (* the x axis of the screen, the b vector, is horizontal and
      orthogonal to the gaze *)
   b[1] := +gaze[2];
   b[2] := -gaze[1];
   b[3] := 0;
   db := sqrt(b[1]*b[1] + b[2]*b[2] + b[3]*b[3]);

   (* check for top view case and correct if so: *)
   if db = 0.0 then begin
      db := 1;
      b[1] := 1;
      b[2] := 0;
   (* b[3] := 0; already from above *)
   end
   else for j := 1 to 3 do b[j] := b[j]/db; (* make b a unit vector *)

   (* now that the gaze is a unit vector, and we have constructed
      the x axis in the b vector also as a unit vector, the cross
      product of these two will generate the y axis as a unit
      vector, c: *)
   c[1] := +(b[2]*gaze[3] - gaze[2]*b[3]);
   c[2] := -(b[1]*gaze[3] - gaze[1]*b[3]);
   c[3] := +(b[1]*gaze[2] - gaze[1]*b[2]);

   (* now normalize both b and c vectors to be of size 1/smag *)
   for j := 1 to 3 do begin
      b[j] := b[j]/smag;
      c[j] := c[j]/smag;
   end
end;
(* end module pic.3d.view version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.3d.makescreen *)
procedure makescreen(vx,vy,vz, gx,gy,gz, smagnitude: real; var s: screen);
(* create the screen s based on the viewing location (vx,vy,vz)
and the direction of gaze (gz,gy,gz).  The screen size is scaled by
smagnitude; doubling smagnitude should double the size of the scene. *)

(* This routine makes creation of the screen very simple for the user.
One need not look at the view routine. *)
begin
   s.v[1] := vx;
   s.v[2] := vy;
   s.v[3] := vz;
   s.g[1] := gx;
   s.g[2] := gy;
   s.g[3] := gz;
   with s do view(v,g,smagnitude, a,b,c);
   s.smag := smagnitude;
   s.range := 1/smagnitude
end;
(* end module pic.3d.makescreen version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.3d.project3d *)
procedure project3d(x,y,z: real; s: screen; var xscreen,yscreen: real);
(* project the point (x,y,z) onto the screen s, to find the screen
coordinates (xscreen and yscreen). *)

(* This routine simplifies the projection function for the user. *)
var
   o: threevector; (* for passing the values to d32 *)
begin
   o[1] := x;
   o[2] := y;
   o[3] := z;
   with s do d32(o,a,b,c,v,xscreen,yscreen);
end;
(* end module pic.3d.project3d version = 2.63; (@ of dops.p 1991 November 2 *)
(* ********************************************************************** *)
(* end module pic.3d.package version = 2.63; (@ of dops.p 1991 November 2 *)


(* ********************************************************************** *)
(* ********************************************************************** *)
(* ********************************************************************** *)

(* begin module pic.3d.test.fun *)
function fun(r: real): real;
(* a function to plot *)
begin
   fun := 3/(1+r*r/2)
end;
(* end module pic.3d.test.fun version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module pic.3d.test.test3d *)
procedure test3d(var afile: text);
(* test three dimensional graphics *)
var
   s: screen; (* the screen on which to project the 3d image *)

   xscreen, yscreen: real; (* location on the screen corresponding to the
      projection of o onto the screen defined by v,a,b,c *)
   oldxscreen,oldyscreen: real; (*  the previous valuse of xscreen and yscreen *)

   (* definition of a spiral *)
   dr: real; (* change in r *)
   dtheta: real; (* change in theta *)
   r: real; (* radius of the current position *)
   radius: real; (* the radius of the spiral *)
   theta: real; (* angle of the current position *)
   thickness: real; (* spacing between spiral arms *)
   steps: real; (* number of steps around a circle of the spiral *)

   x,y,z: real; (* the location in three space *)
begin

   makescreen(5.0,5.0,5.0, -1.0,-1.0,-1.0, 5.0, s);
   
   r := 0;
   theta := 0;
   steps := 15;
   thickness := 0.1;
   radius := 2.0;
   dr := thickness/steps;
   dtheta := 2 * pi / steps;

   x := 0;
   y := 0;
   z := fun(r);
   project3d(x,y,z, s, oldxscreen,oldyscreen);
   mover(afile,oldxscreen,oldyscreen); (* premove to the startpoint of the graph *)

   while r < radius do begin
      r := r + dr;
      theta := theta + dtheta;
      polrec(r,theta,x,y);

      z := fun(r);
      project3d(x,y,z, s, xscreen,yscreen);

      (* draw a line from where we where to the new place *)
      liner(afile, xscreen - oldxscreen, yscreen - oldyscreen);
      oldxscreen := xscreen;
      oldyscreen := yscreen;
   end;
end;
(* end module pic.3d.test.test3d version = 2.63; (@ of dops.p 1991 November 2 *)

(* begin module da3d.getfromto *)
procedure getfromto(var da: text; var f,t: integer);
(* get the from (f) and to (t) range of the da file *)
var c: char; (* for skipping the * on the line *)
begin
   reset(da);
   readln(da); (* skip id line *)
   readln(da); (* skip book line *)
   readln(da,c, f,t); (* skip book line *)
writeln(output,f:1,' ',t:1);
end;
(* end module da3d.getfromto *)

(* begin module da3d.readdata *)
procedure readdata(var da: text; datacolumn: integer;
          var fromrange, torange: integer;
          var data: datagrid);
(* read x,y (columns 1 and 2) and z data from column datacolumn of file da. *)
var
   idcharacter: char; (* identification character for reading data *)
   skipc: char; (* for skipping characters *)
   column: integer; (* for reading through columns *)
   x,y: integer; (* location on the base plane *)
   z: real; (* the location in three space *)
begin
   (* obtain the from-to range *)
   getfromto(da,fromrange,torange);

   if torange - fromrange > maxgrid - mingrid then begin
      writeln(output,'range (from ',fromrange:1,
                     ' to ',torange:1,') is bigger than allowed',
                     ' by constant maxgrid (',maxgrid:1,')');
      halt;
   end;

   (* clear grid *)
   for x := mingrid to maxgrid do begin
      for y := mingrid to maxgrid do begin
         data[x,y] := -maxint;
      end;
   end;

   while not eof(da) do begin
      read(da,idcharacter);
      if idcharacter = 'i' then begin
         read(da,skipc,skipc,skipc); (* skip characters *)
         read(da,x,y);
         (* skip numbers until the z is read: *)
         for column := 5 to datacolumn do read(da, z);
         readln(da);

         if (x < mingrid) or (x > maxgrid) or
            (y < mingrid) or (y > maxgrid) then begin
            writeln(output,'point (',x:1,',',y:1,') is out of range',
                           mingrid:1,' to ',maxgrid:1);
            halt;
         end;

{zzz}
         data[x,y] := z;
      end
      else readln(da);
   end
end;
(* end module da3d.readdata *)

(* begin module da3d.grid *)
procedure grid(var afile: text;
               data: datagrid;
               dodiagonal: boolean;
               fromrange, torange: real;
               xlocation: real;
               ylocation: real;
               zlocation: real;
               magnify: real;
               xmagnify: real;
               ymagnify: real;
               zmagnify: real
               );
(* create three dimensional graphics from the data *)
var
   s: screen; (* the screen on which to project the 3d image *)

   xscreen, yscreen: real; (* location on the screen corresponding to the
      projection of o onto the screen defined by v,a,b,c *)
   oldxscreen,oldyscreen: real; (*  the previous valuse of xscreen and yscreen *)

   x,y: integer; (* location on the base plane *)
   z: real; (* the location in three space *)
   xold,yold: real; (* the old location *)

   changed: boolean; (* if true, the data has changed suddenly;
      use a move to the current coordinate rather than drawing a line *)

   xhold,yhold: real; (* screen coordinates for temporary markings *)

procedure drawit;
(* draw to the point x,y,z *)
begin
         z := data[x,y];
         if (z > -maxint) then if (dodiagonal and (x = y)) or (x <> y)
         then begin
            writeln(afile,'% (',x:3,', ',y:3,
                  ', ',z:10:5,')');

            (* any jump in coordinates restarts the drawing line *)
            if (abs (x - xold) > 1) or
               (abs (y - yold) > 1)
            then changed := true
            else changed := false;

            if changed then begin (* label axis *)
               project3d(x*xmagnify,y*ymagnify,0.0, s, xhold,yhold);
               movea(afile, xhold, yhold);
               project3d(x*xmagnify,y*ymagnify,-1.0, s, xhold,yhold);
               linea(afile, xhold, yhold);
               writeln(afile, 'gsave s (',x:1,') x grestore');
            end;

            project3d(x*xmagnify,y*ymagnify,z*zmagnify, s, xscreen,yscreen);

            if changed
            then movea(afile, xscreen, yscreen)
            else liner(afile, xscreen - oldxscreen, yscreen - oldyscreen);

            oldxscreen := xscreen;
            oldyscreen := yscreen;
            xold := x;
            yold := y;
            changed := false;
         end
end;

begin

   makescreen( xlocation, ylocation, zlocation,
              -xlocation,-ylocation,-zlocation,
               magnify, s);

   (* outline the area *)
   project3d(fromrange*xmagnify,fromrange*ymagnify,0, s, xscreen,yscreen);
   movea(afile, xscreen, yscreen);

   project3d(fromrange*xmagnify,torange*ymagnify,0, s, xscreen,yscreen);
   linea(afile, xscreen, yscreen);

   project3d(fromrange*xmagnify,torange*ymagnify,2.0, s, xhold,yhold);
   linea(afile, xhold, yhold);
   writeln(afile, 'gsave s (',2/(zmagnify*magnify):4:2,' bits) x grestore');
   linea(afile, xscreen, yscreen);

   project3d(torange*xmagnify,torange*ymagnify,0, s, xscreen,yscreen);
   linea(afile, xscreen, yscreen);

   project3d(fromrange*xmagnify,fromrange*ymagnify,0, s, xscreen,yscreen);
   linea(afile, xscreen, yscreen);

   xold := -maxint;
   yold := -maxint;

   for x := mingrid to maxgrid do begin
      for y := mingrid to maxgrid do begin
         drawit;
      end
   end;

   changed := true;

   for y := mingrid to maxgrid do begin
      for x := mingrid to maxgrid do begin
         drawit;
      end
   end;
   
end;
(* end module da3d.grid *)
   
(******************************************************************************)

(* begin module da3d.themain *)
procedure themain(var da, da3dp, scene: text);
(* the main procedure of the program *)
var
   data: datagrid; (* store z values in x,y grid *)
   dd: char; (* character for reading value for setting dodiagonal *)
   dodiagonal: boolean; (* do the diagonal in the graph *)
   fromrange, torange: integer; (* range of the data *)
   magnify: real; (* factor by which to magnify the scene *)
   xmagnify: real; (* factor by which to magnify zscale *)
   ymagnify: real; (* factor by which to magnify zscale *)
   zmagnify: real; (* factor by which to magnify zscale *)
   xlocation: real; (* x location *)
   ylocation: real; (* y location *)
   zlocation: real; (* z location *)
   datacolumn: integer; (* column of da to use *)
   horizontal: real; (* zero coordinate of screen in cm *)
   vertical: real; (* zero coordinate of screen in cm *)
   cmperinch: real; (* conversion factor from cm to inches *)

procedure readda3d(var f: text);
begin
      (* obtain parameters *)
      readln(f,horizontal);
      readln(f,vertical);
      readln(f,xlocation);
      readln(f,ylocation);
      readln(f,zlocation);
      readln(f,magnify);
      readln(f,xmagnify);
      readln(f,ymagnify);
      readln(f,zmagnify);
      readln(f,datacolumn);
      if (datacolumn <= 4) or (datacolumn > 15) then begin
         writeln(output,'data column must be >4 and < 16');
         halt
      end;
      readln(f,dd);
      dodiagonal := (dd = 't');
end;

procedure writeda3d(var f: text);
begin
   (* obtain parameters *)
   writeln(f,horizontal:picfield:picwidth,' horizontal adjust (cm)');
   writeln(f,vertical:picfield:picwidth,' vertical adjust (cm)');
   writeln(f,xlocation:picfield,' xlocation (bases)');
   writeln(f,ylocation:picfield,' ylocation (bases)');
   writeln(f,zlocation:picfield,' zlocation (bases)');
   writeln(f,magnify:picfield:picwidth,' magnify');
   writeln(f,xmagnify:picfield:picwidth,' xmagnify');
   writeln(f,ymagnify:picfield:picwidth,' ymagnify');
   writeln(f,zmagnify:picfield:picwidth,' zmagnify');
   writeln(f,datacolumn:picfield,' data column');
   if (datacolumn <= 4) or (datacolumn > 15) then begin
      writeln(output,'data column must be >4 and < 16');
      halt
   end;
   writeln(f,dd, ' dodiagonal');
end;

begin
   writeln(output,'da3d ',version:4:2);

   rewrite(scene);
   writeln(scene,'%! da3d ',version:4:2);
   writeln(scene);

   reset(da3dp);
   if eof(da3dp)
   then begin
      test3d(scene); (* test routine *)
      startpic(scene,81.0,3.0,5.0,'c');
   end
   else begin
      readda3d(da3dp);
      writeda3d(output);

      cmperinch := 2.54;
      startpic(scene,72.0,horizontal/cmperinch,vertical/cmperinch,'c');

      readdata(da, datacolumn, fromrange, torange, data);

      grid(scene,
           data, dodiagonal,
           fromrange, torange,
           xlocation, ylocation, zlocation,
           magnify,
           xmagnify, ymagnify, zmagnify
           );

   end;

   stoppic(scene);
end;
(* end module da3d.themain *)

begin
   themain(da,da3dp,scene);
1: end.