%!PS

% Exploring cubic catenary fits
% ==================================
% by  Don Lancaster

%  Copyright c 2006 by Don Lancaster & Synergetics, Box 809, Thatcher, AZ, 85552
%  (928) 428-4073 Email: don@tinaja.com Website: http://www.tinaja.com
%  Consulting services available http://www.tinaja.com/info01.html

%  All commercial rights and all electronic media rights ~fully~ reserved.
%  Linking usually welcome. Reposting expressly forbidden. Version 1.1 

% preliminary catenary exploration


%  An example of using my Gonzo Utilities to do a custom engineering graph.
%  Detailed example does a scatterplot on a 3x3 log graph.

%  IMPORTANT NOTE: Don Lancaster's file gonzo.ps is required for this program. 
%  This may be freely downloaded from http://www.tinaja.com/post01.asp#gonzo
%  To activate after obvious location mods, uncomment ONE of the following two lines:

   (C:\\Documents and Settings\\don\\Desktop\\gonzo\\gonzo.ps) run  % use internal gonzo
%  (A:\\gonzo.ps) run  % use external gonzo

%  Plus the next line when Gonzo is to be activated...
   gonzo begin ps.util.1 begin nuisance begin

%  NOTE THAT ALL PS FILENAME STRINGS !!!DEMAND!!! DOUBLE REVERSE SLASHES.

%  To use this program, you modify your values and save it as an ordinary ASCII
%  textfile. That file is then sent to Acrobat Distiller for coversion to a .PDF
%  file format.

mark


% Find y = cosh (x/a)  given a and x on stack ....
% y = a/2 [e^(x/a) + e^(-x/a)]


/ee 2.718281828 store

/cosh { /x1 exch store
        /a1 exch store
        x1 a1 div dup
       ee exch exp       % find e to the xponent
       exch neg
       ee exch exp
       add
       a1 mul
       2 div } store
       






/makeacat {
0 0 mt

-1 0.01 1 {/ii exch store

aa ii cosh

ii 1 add 10 mul exch
10 mul lineto 

} for

1 0 1 setrgbcolor
line1 stroke } store







/makeacat1 {
0 0 mt

-3 0.01 3 {/ii exch store

aa ii cosh

ii 3 add 3.3333 mul exch
5 mul lineto 

} for

1 0 1 setrgbcolor
line3 stroke } store



%%%%%%%%%%%%%%% begin simple plot a = 1

50 50 10 setgrid
20 20 showgrid


0 0 mt 20 pu 20 pr 20 pd closepath 
0.7 1 1 setrgbcolor fill

0 1 1 setrgbcolor line2

1 setlinecap 1 setlinejoin

[ {0 0 mt 20 u} 3.33333 7 ] xrpt
[ {0 0 mt 20 r} 2.5 9 ] yrpt

line2


0 0.8 0.8 setrgbcolor

0 -0.4 mt 20.4 pu 20 pr 20.4 pd

-0.4 0 mt 20.4 r
-0.4 20 mt 20.4 r
10 0 mt 0.4 d
0 10 mt 0.4 l

10 0 mt 20 u
0 10 mt 20 r


gsave
0.3 0 mt 20 pu 20 pr 20 pd closepath clip

newpath

/aa 2 store makeacat1
/aa 1 store makeacat1
/aa 0.5 store makeacat1





% try a series fit for a = 1 to check math

%  y = 1 + xx/2 + xxxx/24 + xxxxxx/720 + ...

% 

10 5 mt


black line1 
0 0.01 3 { /xx exch store


1 

xx dup mul 2 div add   % square term

xx xx mul dup mul 24 div add    % fourth term

xx dup mul dup mul   % x4
xx dup mul mul       % x6
720 div   add

5 mul 

xx 3.333333 mul 10 add  exch lineto 


} for

stroke






% same fit at a=2

% 
%       y = 2 + xx/4 + xxxx/192 + xxxxxx/23040


10 10 mt


black line1 
0 0.01 3 { /xx exch store


2
xx dup mul 4 div add   % square term

xx xx mul dup mul 192 div add    % fourth term

xx dup mul dup mul   % x4
xx dup mul mul       % x6
23040 div   add

5 mul 

xx 3.333333 mul 10 add  exch lineto 


} for 

line1 stroke





% same fit at a=.5

%  y = 0.5 + xx  + xxxx/3 + xxxxxx/22.5 + xxxxxxxx/315




10 2.5 mt


black 0 setlinewidth
0 0.01 3 { /xx exch store


0.5
xx dup mul 1 div add   % square term

xx xx mul dup mul 3 div add    % fourth term

xx dup mul dup mul   % x4
xx dup mul mul       % x6
22.5 div   add

% eighth term barely noticable but...

true {xx dup mul dup mul dup mul 
315 div add} if





5 mul 

xx 3.333333333 mul 10 add  exch lineto 


} for 

0 setlinewidth stroke


%%%% attempt to create a slope controlled cubic


/ee 2.718281828 store

/cosh { /x1 exch store
        /a1 exch store
        x1 a1 div dup
       ee exch exp       % find e to the xponent
       exch neg
       ee exch exp
       add
       a1 mul
       2 div } store


%  dy/dx =  2x  + x^3 + 0.222222 x^5 + .022222 x^7  at a = 0.5 

/findslopea.5 { /xx exch store

  xx 2 mul
  xx dup mul xx mul add

  xx dup mul xx mul  xx dup mul add  0.222222222 mul add
  xx dup mul xx mul dup mul xx mul 0.02222222 mul add

1 atan

} store

% above was WRONG


/findslopea.5 { /xx exch store     % hopefully corrected

  xx 2 mul
  xx dup mul xx mul    4 mul 3 div  add

  xx dup mul xx mul  xx dup mul add     4 mul 15 div   add
  xx dup mul xx mul dup mul xx mul    8 mul 315 div     add

1 atan

} store


% try generalized slope finder. needs aa defined


%  slope = 2x/2!a + 4xxx/4!aaa + 6xxxxx/6!aaaaa + 8xxxxxxx/8!aaaaaaa + ignore
%  slope = x/1!a + xxx/3!aaa + xxxxx/5!aaaaa + xxxxxxx/7!a

%  slope = x/a + xxx/6aaa + xxxxx/120aaaaa + xxxxxxx/5040aaaaaaa

/findslope { /xx exch store      

     xx aa div

     xx dup mul xx mul     aa dup mul aa mul  div  6 div        add

     xx dup mul dup mul xx mul  
     aa dup mul dup mul aa mul  div                120 div      add
 
     xx dup mul dup mul  
     xx dup mul xx mul mul    

     aa dup mul dup mul
     aa dup mul aa mul mul     div                 5040 div     add

     1 atan

(\n\n used new findslope.\n\n) print flush

          } store

         
% calculation is wildly wrong



/findslopea.5 {findslope}  store  % try generic



/fitcat.5 { /extent exch store      % exit angle enthuasiasm
           /intent exch store        % entry angle enthuasiasm
           /x3 exch store            % exit x value
           /x0 exch store           % entry x value

           /y0 0.5 x0 cosh store    % entry y value
           /y3 0.5 x3 cosh store    % exit y value

           /entang x0                % find entry angle 
           findslopea.5  store           

           /x1 x0 entang cos        % find first control points
               intent mul
               add store
           /y1 y0 entang sin
               intent mul
               add store 

        

           /extang x3                % find exitangle 
           findslopea.5  store


(\n\n findslopa.5 on exit is now ) extang
20 string cvs mergestr (.\n\n) mergestr print flush



 


           /x2 x3 extang cos        % find first control points
               extent mul neg
               add store
           /y2 y3 extang sin
               extent mul neg 
               add store  


% save normalized






  /x00 x0 3.33333 mul 10 add store
  /x01 x1 3.33333 mul 10 add store
  /x02 x2 3.33333 mul 10 add store
  /x03 x3 3.33333 mul 10 add store

  /y00 y0 5 mul store
   /y01 y1 5 mul store
   /y02 y2 5 mul store
   /y03 y3 5 mul store




        true { x00 y00 mt dot    %  show dots
         x01 y01 mt dot
         x02 y02 mt dot
         x03 y03 mt dot} if
         

         x00 y00 mt
         x01 y01 x02 y02 x03 y03 curveto


          } store





0 0 1 setrgbcolor 0 setlinewidth




% 0 1 0.6 0.4  fitcat.5  good fit!

% 0 1.386  1.10 0.10  fitcat.5      % 1.10 0.10 is ok


/aa 0.5 store  % reset a

0 1.386     0.837 1.64 fitcat.5

0 setlinewidth stroke

grestore % remove clip


0.2 0.2 0.2 setrgbcolor % ease off on black harsh

/font1 /Helvetica 1 gonzofont
/kern 0 store

font1

0 -1.7 (X|j=|j-1) cc
10 -1.7 (X|j=|j0) cc
20 -1.7 (X|j=|j1) cc

-1 -0.3 (Y|j=|j0) cr
-1  9.7 (Y|j=|j2) cr
-1 19.7 (Y|j=|j4) cr

10 10.6 (a|j=|j2) cc
10 5.6 (a|j=|j1) cc


/kern 0.15 store

10 1.2 (a= 0|k.|k5) cc

showpage


% EOF

% debugging - not using slope info


% error analysis   we used the wrong slope at the end. Any chance it
% matches the right one?

1.386 findslopea.5 ==   % returns 81.4752 degrees  


% old  dy/dx =  2x  +      x^3 + 0.222222 x^5 + .022222 x^7  at a = 0.5 
% new  dy/dx =  2x  + (4/3)x^3 +   (4/15) x^5 + (8/315)  x^7  at a = 0.5 



/betterfindslopea.5 { /xx exch store

  xx 2 mul
  xx dup mul xx mul    4 mul 3 div  add

  xx dup mul xx mul  xx dup mul add     4 mul 15 div   add
  xx dup mul xx mul dup mul xx mul    8 mul 315 div     add

1 atan

} store


1.386 betterfindslopea.5 ==    


%%%%%%%%%%%%% try to plot error dist


/findAH { 
/A x3 x2 3 mul sub x1 3 mul add x0 sub store
/E y3 y2 3 mul sub y1 3 mul add y0 sub store
/B x2 3 mul x1 6 mul sub x0 3 mul add  store
/F y2 3 mul y1 6 mul sub y0 3 mul add  store
/C x1 3 mul x0 3 mul sub store
/G y1 3 mul y0 3 mul sub store
/D x0 store
/H y0 store
        } store


findAH   % assumes we still have valid points

A ==
B ==
C ==
D ==




/showdiff {    /rmserr 0 store

0 .01 1.001 {/tt exch store

         D  
         C     tt                   mul add
         B     tt dup mul           mul add
         A     tt dup mul tt mul    mul add  
         /xxx exch store

         H
         G tt mul add
         F tt dup mul mul add
         E tt dup mul tt mul mul add 
         /yyy exch store

  0.5 xxx cosh

 /wanty exch store


% (\n\n) print flush
 %  xxx ==
 %  yyy ==
  % wanty ==
  wanty yyy sub

 dup /curerr exch store

 ==

/rmserr rmserr curerr dup mul add store




   } for


(\n\nrms error is ) rmserr sqrt 100 div 20 string cvs mergestr (.\n\n) mergestr print



   } store




showdiff


(\n\n next pass             sdgsdfg\n\n) print flush

0 1.386     0.837 .004 add           1.64   0.002 add               fitcat.5

findAH

showdiff

% .024 add    rms error is 0.00149825.

% .014 add rms error is    0.000758119.
% .012 add rms error is    0.000616709.

% .010 add                 0.000481645.

% .007 add    rms error is 0.000308943.
% .006 add    rms error is 0.000269744
% .005 add    rms error is 0.000247894.   <------------------------------------  
% .004 add    rms error is 0.000247988                           


% .000 add    rms error is 0.000417987.


% .012 sub    rms error is 0.00125967.
% .013 sub    rms error is 0.00133273.
% .014 sub is rms error is 0.00140589.  <---
% .015 sub rms error is    0.00147912.


% now adjust high end

% 0.000 add   rms error is 0.000247894.
% 0.001 add   rms error is 0.000246958
% 0.002 add   rms error is 0.000246675
% 0.003 add   rms error is 0.000247216

% 0.006 add   rms error is 0.000253526.

% best fit appears to be 0.841 and 1.642

% not clear whether local or global miminum.