program coscurve(coscurvep, output);
(* coscurve: compute the length of the cosine curve

  Dr. Thomas D. Schneider
  National Cancer Institute
  Laboratory of Experimental and Computational Biology
  Frederick, Maryland  21702-1201
  toms@ncifcrf.gov
  permanent email: toms@alum.mit.edu
  http://www.lecb.ncifcrf.gov/~toms/

  National Cancer Institute
  Laboratory of Experimental and Computational Biology

 *)

label 1; (* end of program *)

const
(* begin module version *)
version = 1.05; (* of coscurve.p 2000 Oct 31
2000 Oct 31, 1.00: origin
*)
updateversion = 1.03; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.coscurve *)
(*
name
   coscurve: compute the length of the cosine curve

synopsis
   coscurve(coscurvep: in, output: out)

files

   coscurvep:  parameters to control the program.  The file must contain the
      following parameters, one per line:

      parameterversion: The version number of the program.  This allows the
         user to be warned if an old parameter file is used.

      second line: intervals (integer):  number of intervals to use

      third line: wid (integer):  width of the output real numbers (in characters)

      fourth line:  printcontrol: the first character determines how or
         whether to print the values during the summation.

         P: print x, y, dy, L, L/2pi at each step
         p: print              L/2pi at each step
         m: print              L/2pi at each successively larger power of 10
         -: no printing (fastest compute)

   output: messages to the user

description

   The cosine wave on a sequence logo can have dashes.  These are implemented
   in PostScript using the PostScript dash function.  To make the dashes line
   up with the wavelength the length of the curve must be known.  So the
   problem is to compute length of y = cos(x) from 0 to 2 pi.  Setting ds =
   sqrt(dx^2 +dy^2) (Thomas, pages 203-204) we find that the total length is
   L = integral sqrt(1 + sin^2(x)) dx.  This is not solvable in terms of
   elementary functions (Thomas, page 308) and so must be determined
   numerically.  That is the purpose of this program.

   Rather than integrate the function, it is easier just to sum the ds for
   the interval 0 to 2 pi.  The value L/(2pi) is reported so that this need
   not be computed in the final program, which uses the constant.

   The user may then multiply L/(2pi) by the amplitude and wavelength of a
   given cosine or sine wave to get its contour length.

examples

for cocurvep being:

1.03      version of coscurve that this parameter file is designed for.
100000000 number of intervals to use
20        width of output real numbers in characters
m 2       printing: P: all, p: just L/2pi, m: multiples of k, - no printing.

the output is:

coscurve 1.05
pi =  3.14159265358979312
intervals = 100000000
                   2, L/(2pi) =  1.18544706105728359
                   4, L/(2pi) =  1.18544706105728359
                   8, L/(2pi) =  1.20642353304458783
                  16, L/(2pi) =  1.21356557113189623
                  32, L/(2pi) =  1.21539466952667086
                  64, L/(2pi) =  1.21585359980789032
                 128, L/(2pi) =  1.21596843563495249
                 256, L/(2pi) =  1.21599715104708350
                 512, L/(2pi) =  1.21600433030360899
                1024, L/(2pi) =  1.21600612514295325
                2048, L/(2pi) =  1.21600657385437816
                4096, L/(2pi) =  1.21600668603233197
                8192, L/(2pi) =  1.21600671407678718
               16384, L/(2pi) =  1.21600672108796770
               32768, L/(2pi) =  1.21600672284060907
               65536, L/(2pi) =  1.21600672327877302
              131072, L/(2pi) =  1.21600672338823634
              262144, L/(2pi) =  1.21600672341558114
              524288, L/(2pi) =  1.21600672342506755
             1048576, L/(2pi) =  1.21600672342277272
             2097152, L/(2pi) =  1.21600672342234084
             4194304, L/(2pi) =  1.21600672341971094
             8388608, L/(2pi) =  1.21600672346093108
            16777216, L/(2pi) =  1.21600672339648463
            33554432, L/(2pi) =  1.21600672338265259
            67108864, L/(2pi) =  1.21600672333825899
L =  7.64039557751054144
L/(2pi) =  1.21600672333825899
==========

documentation

   George B. Thomas, Jr.
   Calculus and Analytic Genometry, 4th Ed
   Addison-Wesley Publishing Co.
   Reading, MA. 1968.

see also

   {Parameter file: } coscurvep
   {The program that uses this number:} makelogo.p

author

   Thomas Dana Schneider

bugs

technical notes

*)
(* end module describe.coscurve *)

var coscurvep: text; (* parameter file *)

(* 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 = 'delmod 6.16 84 mar 12 tds/gds'; *)

(* begin module coscurve.themain *)
procedure themain(var coscurvep: text);
(* the main procedure of the program *)
var
   dec: integer; (* decimals of output numbers *)
   dx: real; (* change in x *)
   dy: real; (* change in y *)
   intervals: integer; (* number of intervals to compute over *)
   L: real; (* total length *)
   multiple: integer; (* multiple increase density of points by *)
   parameterversion: real; (* parameter version number *)
   pi: real; (* 3.1415... *)
   printcontrol: char; (* control of printing *)
   printing: boolean; (* if true, print intermediate values *)
   twopi: real; (* twice pi *)
   target: integer; (* next target number of intervals at which to print at *)
   wid: integer; (* width of output numbers *)

procedure compute(intervals: integer);
(* compute the length for a given number of intervals *)
var
   i: integer; (* index to intervals *)
   x: real; (* x axis value *)
   y: real; (* y axis value *)
   xold: real; (* previous value of x *)
   yold: real; (* previous value of y *)
begin
   dx := twopi/intervals;
   if printcontrol in ['P','p']
   then writeln(output, 'dx = ', dx:wid);
   L := 0.0;
   x := 0.0;
   xold := 0.0;
   yold := cos(xold);
   for i := 1 to intervals do begin
      x := x + dx;
      y := cos(x);
      dy := y - yold;
      L := L + sqrt(dx*dx + dy*dy);
      if printing then begin
         if printcontrol in ['P','p'] then begin
            if printcontrol = 'P' then begin
               write(output, '  x = ',  x:wid:dec);
               write(output, '  y = ',  y:wid:dec);
               write(output, ' dy = ', dy:wid:dec);
               write(output, '  L = ',  L:wid:dec);
               writeln(output);
            end;
            writeln(output, '  L/(2pi) = ',  L/twopi:wid:dec);
         end;
      end;
      xold := x;
      yold := y;
   end;
end;

begin
   writeln(output,'coscurve ',version:4:2);
   reset(coscurvep);
   readln(coscurvep, parameterversion);
   if round(100*parameterversion) < round(100*updateversion) then begin
      writeln(output, 'You have an old parameter file!');
      halt
   end;

   readln(coscurvep, intervals);
   readln(coscurvep, wid);

   read(coscurvep, printcontrol);
   if printcontrol = 'm'
   then begin
      read(coscurvep, multiple);
      if multiple <= 1 then begin
         writeln(output, 'multiple must be >= 1 ');
         halt;
      end;
   end
   else multiple := 1;
   readln(coscurvep);

   if not (printcontrol in ['P','p','m','-']) then begin
      writeln(output, 'printcontrol must be one of: Ppm-');
      halt;
   end;

   printing := (printcontrol <> '-');
   dec := wid - 3;

   pi := 4.0 * arctan(1.0);
   twopi := 2.0 * pi;
   writeln(output,'pi = ',pi:wid:dec);

   writeln(output, 'intervals = ', intervals:1);

   if printcontrol = 'm'
   then begin
      target := multiple;
      while target <= intervals do begin
         compute(target);
         writeln(output, target:wid, ', L/(2pi) = ',  L/twopi:wid:dec);
         target := multiple*target;
      end;
   end
   else begin
      compute(intervals);
   end;

   writeln(output, 'L = ',L:wid:dec);
   writeln(output, 'L/(2pi) = ',  L/twopi:wid:dec);
   writeln(output, '==========');

end;
(* end module coscurve.themain *)

begin
   themain(coscurvep);
1: end.