home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 509.lha / ATS_v2.0 / ats2.cp < prev    next >
Text File  |  1991-05-06  |  5KB  |  121 lines
  1. //        ret.terms[i]=complex(0.0,0.0);
  2. //
  3. //    C++ Routines to manipulate Taylor series for solving a differential
  4. //  equation.
  5. //
  6. //  By Dennis J. Darland 
  7. //  4/3/91
  8. //  Public Domain
  9. //  Version 1.6
  10. //  Many enhancements expected.
  11. //  PROVIDED WITH NO WARRANTY.
  12. //  
  13. #include "ats.h"
  14. complex radius_cp(t_s & y,complex & order)
  15.     {
  16.     register int m;
  17.     complex rm0,rm1,rm2,rm3,rm4,rcs,c;
  18.     complex t1,t2,t3,t4;
  19.         
  20.     m = y.last_term - 1;
  21.     rm0 = y.terms[m]/y.terms[m-1];
  22.     cout << "rm0 = " << rm0 << "\n";
  23.     rm1 = y.terms[m-1]/y.terms[m-2];
  24.     cout << "rm1 = " << rm1 << "\n";
  25.     rm2 = y.terms[m-2]/y.terms[m-3];
  26.     cout << "rm2 = " << rm2 << "\n";
  27.     rm3 = y.terms[m-3]/y.terms[m-4];
  28.     cout << "rm3 = " << rm3 << "\n";
  29.     rm4 = y.terms[m-4]/y.terms[m-5];
  30.     cout << "rm4 = " << rm4 << "\n";
  31.     t1 = complex(12.0,0.0)*rm2*rm1*rm1*rm4;
  32.     t1 = t1 -complex(6.0,0.0)*rm3*rm3*rm2*rm4;
  33.     t1 = t1 -complex(10.0,0.0)*rm1*rm1*rm2*rm3;
  34.     t1 = t1 +complex(8.0,0.0)*rm2*rm2*rm3*rm4;
  35.     t1 = t1 +complex(10.0,0.0)*rm2*rm2*rm1*rm3;
  36.     t1 = t1 +complex(18.0,0.0)*rm1*rm3*rm3*rm4;
  37.     t1 = t1 +complex(25.0,0.0)*complex(double(m),0.0)*rm1*rm3*rm2*rm4;
  38.     t1 = t1 -complex(5.0,0.0)*complex(double(m),0.0)*rm0*rm2*rm1*rm3;
  39.     t1 = t1 +complex(12.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm3;
  40.     t1 = t1 -complex(7.0,0.0)*complex(double(m),0.0)*rm2*rm2*rm1*rm3;
  41.     t1 = t1 -complex(32.0,0.0)*rm2*rm1*rm3*rm4;
  42.     t1 = t1 -complex(8.0,0.0)*complex(double(m),0.0)*rm2*rm2*rm3*rm4;
  43.     t1 = t1 +complex(5.0,0.0)*complex(double(m),0.0)*rm3*rm3*rm2*rm4;
  44.     t1 = t1 -complex(12.0,0.0)*complex(double(m),0.0)*rm3*rm3*rm4*rm1;
  45.     t1 = t1 -complex(15.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm4;
  46.     t1 = t1 +complex(8.0,0.0)*complex(double(m),0.0)*rm0*rm2*rm1*rm4;
  47.     t1 = t1 -complex(3.0,0.0)*complex(double(m),0.0)*rm0*rm1*rm3*rm4;
  48.     t1 = t1 +rm1*rm3*rm4*rm0*complex(double(m),0.0)*complex(double(m),0.0);
  49.     t1 = t1 +rm2*rm1*rm3*rm0*complex(double(m),0.0)*complex(double(m),0.0);
  50.     t1 = t1 -complex(2.0,0.0)*rm2*rm1*rm4*rm0*complex(double(m),0.0)*complex(double(m),0.0);
  51.     t1 = t1 -complex(2.0,0.0)*rm1*rm1*rm2*rm3*complex(double(m),0.0)*complex(double(m),0.0);
  52.     t1 = t1 +complex(3.0,0.0)*rm2*rm1*rm1*rm4*complex(double(m),0.0)*complex(double(m),0.0);
  53.     t1 = t1 +rm2*rm2*rm1*rm3*complex(double(m),0.0)*complex(double(m),0.0);
  54.     t1 = t1 -complex(5.0,0.0)*rm2*rm1*rm3*rm4*complex(double(m),0.0)*complex(double(m),0.0);
  55.     t1 = t1 +complex(2.0,0.0)*rm1*rm3*rm3*rm4*complex(double(m),0.0)*complex(double(m),0.0);
  56.     t1 = t1 +complex(2.0,0.0)*rm2*rm2*rm3*rm4*complex(double(m),0.0)*complex(double(m),0.0);
  57.     t1 = t1 -rm3*rm3*rm2*rm4*complex(double(m),0.0)*complex(double(m),0.0);
  58.     t2 = -complex(5.0,0.0)*complex(double(m),0.0)*rm1*rm3*rm2*rm4;
  59.     t2 = t2 +complex(3.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm4;
  60.     t2 = t2 +complex(10.0,0.0)*rm2*rm1*rm3*rm4;
  61.     t2 = t2 -complex(3.0,0.0)*rm2*rm1*rm1*rm4;
  62.     t2 = t2 +complex(2.0,0.0)*complex(double(m),0.0)*rm2*rm2*rm3*rm4;
  63.     t2 = t2 -complex(4.0,0.0)*rm2*rm2*rm3*rm4;
  64.     t2 = t2 +complex(2.0,0.0)*complex(double(m),0.0)*rm3*rm3*rm4*rm1;
  65.     t2 = t2 -complex(double(m),0.0)*rm3*rm3*rm2*rm4;
  66.     t2 = t2 -complex(6.0,0.0)*rm1*rm3*rm3*rm4;
  67.     t2 = t2 +complex(3.0,0.0)*rm3*rm3*rm2*rm4;
  68.     t2 = t2 -complex(2.0,0.0)*complex(double(m),0.0)*rm0*rm2*rm1*rm4;
  69.     t2 = t2 +complex(double(m),0.0)*rm0*rm1*rm3*rm4;
  70.     t2 = t2 +complex(double(m),0.0)*rm0*rm2*rm1*rm3;
  71.     t2 = t2 -complex(2.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm3;
  72.     t2 = t2 +complex(2.0,0.0)*rm1*rm1*rm2*rm3;
  73.     t2 = t2 +complex(double(m),0.0)*rm2*rm2*rm1*rm3;
  74.     t2 = t2 -complex(2.0,0.0)*rm2*rm2*rm1*rm3;
  75.     order = -t1/(complex(2.0,0.0)*t2);
  76.     t3 = rm3*rm4;
  77.     t3 = t3 -complex(4.0,0.0)*rm2*rm4;
  78.     t3 = t3 +rm2*rm1;
  79.     t3 = t3 -complex(4.0,0.0)*rm1*rm3;
  80.     t3 = t3 +complex(3.0,0.0)*rm1*rm4;
  81.     t3 = t3 +complex(3.0,0.0)*rm2*rm3;
  82.     t4 = -complex(5.0,0.0)*complex(double(m),0.0)*rm1*rm3*rm2*rm4;
  83.     t4 = t4 +complex(3.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm4;
  84.     t4 = t4 +complex(10.0,0.0)*rm2*rm1*rm3*rm4;
  85.     t4 = t4 -complex(3.0,0.0)*rm2*rm1*rm1*rm4;
  86.     t4 = t4 +complex(2.0,0.0)*complex(double(m),0.0)*rm2*rm2*rm3*rm4;
  87.     t4 = t4 -complex(4.0,0.0)*rm2*rm2*rm3*rm4;
  88.     t4 = t4 +complex(2.0,0.0)*complex(double(m),0.0)*rm3*rm3*rm4*rm1;
  89.     t4 = t4 -complex(double(m),0.0)*rm3*rm3*rm2*rm4;
  90.     t4 = t4 -complex(6.0,0.0)*rm1*rm3*rm3*rm4;
  91.     t4 = t4 +complex(3.0,0.0)*rm3*rm3*rm2*rm4;
  92.     t4 = t4 -complex(2.0,0.0)*complex(double(m),0.0)*rm0*rm2*rm1*rm4;
  93.     t4 = t4 +complex(double(m),0.0)*rm0*rm1*rm3*rm4;
  94.     t4 = t4 +complex(double(m),0.0)*rm0*rm2*rm1*rm3;
  95.     t4 = t4 -complex(2.0,0.0)*complex(double(m),0.0)*rm1*rm1*rm2*rm3;
  96.     t4 = t4 +complex(2.0,0.0)*rm1*rm1*rm2*rm3;
  97.     t4 = t4 +complex(double(m),0.0)*rm2*rm2*rm1*rm3;
  98.     t4 = t4 -complex(2.0,0.0)*rm2*rm2*rm1*rm3;
  99.     rcs  = -t3/t4;
  100.     c = my_sqrt(rcs)*y.h;
  101.     return(c);
  102.     }
  103. complex radius_xx(t_s & y,complex & order)
  104.     {
  105.     register int m;
  106.     complex rm0,rm1,hdrc,rcs;
  107.         
  108.     m = y.last_term - 1;
  109.     rm0 = y.terms[m]/y.terms[m-1];
  110.     cout << "rm0 = " << rm0 << "\n";
  111.     rm1 = y.terms[m-1]/y.terms[m-2];
  112.     cout << "rm1 = " << rm1 << "\n";
  113.     hdrc = complex(double(m),0.0)*rm0-complex(double(m-1),0.0)*rm1;
  114.     cout << "hdrc = " << hdrc << "\n";
  115.     rcs = y.h/hdrc;
  116.     cout << "rcs = " << rcs << "\n";
  117.     order = complex(double(m),0.0)*rm0*rcs/y.h-(complex(double(m-1),0.0));
  118.     cout << "order = " << order << "\n";
  119.     return(rcs);
  120.     }    
  121.