home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / modellib / class_def.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-02-13  |  21.8 KB  |  1,134 lines
  1. /*
  2.     This file is automatically written by kwrite_f.
  3. */
  4. /* KAOS DYNAMICAL SYSTEM CLASS = class_demo */
  5.  
  6. #include "model.h"
  7.  
  8.  
  9. int d4hammm_init()
  10. {
  11.     title_label = "D4 Nilpotent Hamiltonian --";
  12.  
  13.     mapping_on = 0;
  14.     inverse_on = 1;
  15.     fderiv_on = 0;
  16.     enable_polar = 1;
  17.     enable_period = 0;
  18.     
  19.     var_dim = 4;
  20.     param_dim = 2;
  21.     func_dim = 3;
  22.     (void) malloc_init();
  23.  
  24.     var_label[0] = "x";
  25.     var_label[1] = "y";
  26.     var_label[2] = "z";
  27.     var_label[3] = "w";
  28.     var_polar_label[0] = "r";
  29.     var_polar_label[1] = "rp";
  30.     var_polar_label[2] = "theta";
  31.     var_polar_label[3] = "thetap";
  32.     param_label[0] = "mu";
  33.     param_label[1] = "delta";
  34.     func_label[0] = "Energy";
  35.     func_label[1] = "AngMom";
  36.     func_label[2] = "t";
  37.  
  38.     param[0] = 2;
  39.     param[1] = 1;
  40.  
  41.     var_i[0] = 0;
  42.     var_i[1] = 0;
  43.     var_i[2] = 0;
  44.     var_i[3] = 0;
  45.     var_polar_i[0] = 0;
  46.     var_polar_i[1] = 0;
  47.     var_polar_i[2] = 0;
  48.     var_polar_i[3] = 0;
  49.  
  50.     param_min[0]= -5; param_max[0]= 5;
  51.     param_min[1]= -5; param_max[1]= 5;
  52.  
  53.     var_min[0]= -5; var_max[0]= 5;
  54.     var_min[1]= -5; var_max[1]= 5;
  55.     var_min[2]= -5; var_max[2]= 5;
  56.     var_min[3]= -5; var_max[3]= 5;
  57.  
  58.     var_polar_min[0]= -5; var_polar_max[0]= 5;
  59.     var_polar_min[1]= -5; var_polar_max[1]= 5;
  60.     var_polar_min[2]= -pi; var_polar_max[2]= pi;
  61.     var_polar_min[3]= -5; var_polar_max[3]= 5;
  62.  
  63.     func_min[0]= -1; func_max[0]= 1;
  64.     func_min[1]= -1; func_max[1]= 1;
  65.     func_min[2]= 0; func_max[2]= 10000;
  66.  
  67.     f_p = d4hammm_f;
  68.     func_p = d4hammm_func;
  69. }
  70. /* D4 Hamitonian */
  71. int d4hammm_f(f,index,x,p,t,dim)
  72. int index,dim;
  73. double f[],x[],p[],t;
  74. {
  75.     double v0sq,v2sq;
  76.  
  77.     /* Hamiltonian vector field: 0: dH/dy 2: dH/dw 1:-dH/dx 3:-dH/dz */
  78.     if(index !=2) {
  79.         f[0] = x[1];
  80.         f[2] = x[3];
  81.     }
  82.     if(index !=1){
  83.         v0sq = x[0] * x[0];
  84.         v2sq = x[2] * x[2];
  85.         f[1] = x[0] * (p[0] - (v0sq+v2sq)) + p[1] * x[0] * v2sq;
  86.         f[3] = x[2] *(p[0] - (v0sq+v2sq))+p[1]*x[2]*v0sq; 
  87.     }
  88.     
  89. }
  90. int d4hammm_func(f,x,p,t,dim)
  91. double f[],x[],p[],t;
  92. int dim;
  93. {
  94.         double v0sq,v1sq,v2sq,v3sq;
  95.  
  96.     /* energy and angular momentum are defined in euclidean
  97.     coordinates */
  98.     v0sq = x[0] * x[0];
  99.     v2sq = x[2] * x[2];
  100.     f[0] = 0.5 * ((x[1] * x[1] + x[3] * x[3]) - p[0]*(v0sq+v2sq) + 0.5 *
  101.     (v0sq+v2sq)*(v0sq+v2sq) - p[1] * v0sq* v2sq);
  102.     f[1] =  x[1] * x[2] - x[0] * x[3];
  103.     f[2] = t;
  104.     
  105. }
  106.  
  107. int d4hampp_init()
  108. {
  109.     title_label = "D4 Nilpotent Hamiltonian ++";
  110.  
  111.     mapping_on = 0;
  112.     inverse_on = 1;
  113.     fderiv_on = 0;
  114.     enable_polar = 1;
  115.     enable_period = 0;
  116.     
  117.     var_dim = 4;
  118.     param_dim = 2;
  119.     func_dim = 3;
  120.     (void) malloc_init();
  121.  
  122.     var_label[0] = "x";
  123.     var_label[1] = "y";
  124.     var_label[2] = "z";
  125.     var_label[3] = "w";
  126.     var_polar_label[0] = "r";
  127.     var_polar_label[1] = "rp";
  128.     var_polar_label[2] = "theta";
  129.     var_polar_label[3] = "thetap";
  130.     param_label[0] = "mu";
  131.     param_label[1] = "delta";
  132.     func_label[0] = "Energy";
  133.     func_label[1] = "AngMom";
  134.     func_label[2] = "t";
  135.  
  136.     param[0] = 2;
  137.     param[1] = 1;
  138.  
  139.     var_i[0] = 0;
  140.     var_i[1] = 0;
  141.     var_i[2] = 0;
  142.     var_i[3] = 0;
  143.     var_polar_i[0] = 0;
  144.     var_polar_i[1] = 0;
  145.     var_polar_i[2] = 0;
  146.     var_polar_i[3] = 0;
  147.  
  148.     param_min[0]= -5; param_max[0]= 5;
  149.     param_min[1]= -5; param_max[1]= 5;
  150.  
  151.     var_min[0]= -5; var_max[0]= 5;
  152.     var_min[1]= -5; var_max[1]= 5;
  153.     var_min[2]= -5; var_max[2]= 5;
  154.     var_min[3]= -5; var_max[3]= 5;
  155.  
  156.     var_polar_min[0]= -5; var_polar_max[0]= 5;
  157.     var_polar_min[1]= -5; var_polar_max[1]= 5;
  158.     var_polar_min[2]= -pi; var_polar_max[2]= pi;
  159.     var_polar_min[3]= -5; var_polar_max[3]= 5;
  160.  
  161.     func_min[0]= -1; func_max[0]= 1;
  162.     func_min[1]= -1; func_max[1]= 1;
  163.     func_min[2]= 0; func_max[2]= 10000;
  164.  
  165.     f_p = d4hampp_f;
  166.     func_p = d4hampp_func;
  167. }
  168. /* D4 Hamitonian */
  169. int d4hampp_f(f,index,x,p,t,dim)
  170. int index,dim;
  171. double f[],x[],p[],t;
  172. {
  173.     double v0sq,v2sq;
  174.  
  175.     /* Hamiltonian vector field: 0: dH/dy 2: dH/dw 1:-dH/dx 3:-dH/dz */
  176.     if(index !=2) {
  177.         f[0] = x[1];
  178.         f[2] = x[3];
  179.     }
  180.     if(index !=1){
  181.         v0sq = x[0] * x[0];
  182.         v2sq = x[2] * x[2];
  183.         f[1] = x[0] * (p[0] + (v0sq+v2sq)) + p[1] * x[0] * v2sq;
  184.         f[3] = x[2] * (p[0] + (v0sq+v2sq)) + p[1] * x[2] * v0sq; 
  185.     }
  186.     
  187. }
  188. int d4hampp_func(f,x,p,t,dim)
  189. double f[],x[],p[],t;
  190. int dim;
  191. {
  192.         double v0sq,v1sq,v2sq,v3sq;
  193.  
  194.     /* energy and angular momentum are defined in euclidean
  195.     coordinates */
  196.     v0sq = x[0] * x[0];
  197.     v2sq = x[2] * x[2];
  198.     f[0] = 0.5 * (( x[1] * x[1] + x[3] * x[3]) - p[0]*(v0sq+v2sq)
  199.         + 0.5 * (v0sq+v2sq)*(v0sq+v2sq) - p[1] * v0sq* v2sq);
  200.     f[1] =  x[1] * x[2] - x[0] * x[3];
  201.     f[2] =  t;
  202. }
  203.  
  204. int d4dissmmp_init()
  205. {
  206.     title_label = "D4 Diss. -- +";
  207.  
  208.     mapping_on = 0;
  209.     inverse_on = 1;
  210.     fderiv_on = 0;
  211.     enable_polar = 1;
  212.     enable_period = 0;
  213.     
  214.     var_dim = 4;
  215.     param_dim = 6;
  216.     func_dim = 2;
  217.     (void) malloc_init();
  218.  
  219.     var_label[0] = "x";
  220.     var_label[1] = "y";
  221.     var_label[2] = "z";
  222.     var_label[3] = "w";
  223.     var_polar_label[0] = "r";
  224.     var_polar_label[1] = "rp";
  225.     var_polar_label[2] = "theta";
  226.     var_polar_label[3] = "thetap";
  227.     param_label[0] = "mu";
  228.     param_label[1] = "delta";
  229.     param_label[2] = "epsilon";
  230.     param_label[3] = "nu";
  231.     param_label[4] = "Axzw";
  232.     param_label[5] = "Ayz2";
  233.     func_label[0] = "Energy";
  234.     func_label[1] = "AngMom";
  235.  
  236.     param[0] = 2;
  237.     param[1] = 0;
  238.     param[2] = 0;
  239.     param[3] = 0;
  240.     param[4] = 0;
  241.     param[5] = 0;
  242.  
  243.     var_i[0] = 0;
  244.     var_i[1] = 0;
  245.     var_i[2] = 0;
  246.     var_i[3] = 0;
  247.     var_polar_i[0] = 0;
  248.     var_polar_i[1] = 0;
  249.     var_polar_i[2] = 0;
  250.     var_polar_i[3] = 0;
  251.  
  252.     param_min[0]= -5; param_max[0]= 5;
  253.     param_min[1]= -5; param_max[1]= 5;
  254.     param_min[2]= -5; param_max[2]= 5;
  255.     param_min[3]= -5; param_max[3]= 5;
  256.     param_min[4]= -5; param_max[4]= 5;
  257.     param_min[5]= -5; param_max[5]= 5;
  258.  
  259.     var_min[0]= -5; var_max[0]= 5;
  260.     var_min[1]= -5; var_max[1]= 5;
  261.     var_min[2]= -5; var_max[2]= 5;
  262.     var_min[3]= -5; var_max[3]= 5;
  263.  
  264.     var_polar_min[0]= -5; var_polar_max[0]= 5;
  265.     var_polar_min[1]= -5; var_polar_max[1]= 5;
  266.     var_polar_min[2]= -pi; var_polar_max[2]= pi;
  267.     var_polar_min[3]= -5; var_polar_max[3]= 5;
  268.  
  269.     f_p = d4dissmmp_f;
  270.     func_p = d4dissmmp_func;
  271. }
  272. int d4dissmmp_f(f,index,x,p,t,dim)
  273. int index,dim;
  274. double f[],x[],p[],t;
  275. {
  276.     double v0sq,v2sq;
  277.     if(index !=2) {
  278.         f[0] = x[1];
  279.         f[2] = x[3];
  280.     }
  281.     if(index !=1){
  282.         v0sq = x[0] * x[0];
  283.         v2sq = x[2] * x[2];
  284.         f[1] = x[0] * (p[0]-(v0sq+v2sq)) + p[1] * x[0] * v2sq
  285.             + p[2] * (-(v0sq+v2sq) * x[1] + p[3] * x[1] + p[4] * x[0]*(x[0]*x[1]+x[2]*x[3]) + p[5] * x[1] * v2sq);
  286.         f[3] = x[2] *(p[0] - (v0sq+v2sq))+p[1]*x[2]*v0sq
  287.             + p[2] * (-(v0sq+v2sq) * x[3] + p[3] * x[3] + p[4] * x[2]*(x[0]*x[1]+x[2]*x[3]) + p[5] * x[3] * v0sq);
  288.     }
  289. }
  290. int d4dissmmp_func(f,x,p,t,dim)
  291. double f[],x[],p[],t;
  292. int dim;
  293. {
  294.         double v0sq,v1sq,v2sq,v3sq;
  295.  
  296.     /* energy and angular momentum are defined in euclidean
  297.     coordinates */
  298.     v0sq = x[0] * x[0];
  299.     v2sq = x[2] * x[2];
  300.     f[0] = 0.5 * ((x[1] * x[1] + x[3] * x[3]) - p[0] * (v0sq+v2sq)
  301.         + 0.5 * (v0sq+v2sq) * (v0sq+v2sq) - p[1] * v0sq * v2sq);
  302.     f[1] =  x[1] * x[2] - x[0] * x[3];
  303.     
  304. }
  305. /* Lorenz system */
  306.  
  307. int lorenz_init()
  308. {
  309.     title_label = "Lorenz system";
  310.  
  311.     mapping_on = 0; 
  312.     inverse_on = 1;
  313.     fderiv_on = 0;
  314.     enable_polar = 0;
  315.     enable_period = 0;
  316.  
  317.     var_dim = 3;
  318.     param_dim = 3;
  319.     func_dim = 3;
  320.  
  321.     (void) malloc_init();
  322.  
  323.     var_label[0] = "x";
  324.     var_label[1] = "y";
  325.     var_label[2] = "z";
  326.     param_label[0] = "sigma";
  327.     param_label[1] = "rho";
  328.     param_label[2] = "beta";
  329.     func_label[0] = "t";
  330.     func_label[1] = "x+y";
  331.     func_label[2] = "x-y";
  332.  
  333.     param[0] = 10;
  334.     param[1] = 28;
  335.     param[2] = 2.6666666666666;
  336.     var_i[0] = 0.1;
  337.     var_i[1] = 0.1;
  338.     var_i[2] = 0.1;
  339.  
  340.     param_min[0]= -20; param_max[0]= 20;
  341.     param_min[1]= 0; param_max[1]= 40;
  342.     param_min[2]= -5; param_max[2]= 5;
  343.     var_min[0]= -30; var_max[0]= 30;
  344.     var_min[1]= -30; var_max[1]= 30;
  345.     var_min[2]= -5; var_max[2]= 50;
  346.     func_min[0]= 0; func_max[0]= 100;
  347.     func_min[1]= -60; func_max[1]= 60;
  348.     func_min[2]= -60; func_max[2]= 60;
  349.  
  350.     f_p = lorenz_f;
  351.     func_p = lorenz_func;
  352. }
  353.  
  354. int lorenz_f(f,index,x,p,t,dim)
  355. int index,dim;
  356. double f[],x[],p[],t;
  357. {
  358.     f[0] = p[0] * ( x[1] - x[0] );
  359.     f[1] = p[1] * x[0] - x[1] - x[0] * x[2];
  360.     f[2] = - p[2] * x[2] + x[0] * x[1];
  361.     
  362. }
  363. int lorenz_func(f,x,p,t,dim)
  364. double f[],x[],p[],t;
  365. int dim;
  366. {
  367.     f[0] = t;
  368.     f[1] = x[0] + x[1];
  369.     f[2] = x[0] - x[1];
  370. }
  371. int nlmathieu_init()
  372. {
  373.     title_label = "Nonlinear Mathieu Eq";
  374.  
  375.     mapping_on = 0;
  376.     inverse_on = 1;
  377.     fderiv_on = 0;
  378.     enable_polar = 0;
  379.     enable_period = 0;
  380.  
  381.     var_dim = 3;
  382.     param_dim = 3;
  383.     func_dim = 2;
  384.  
  385.     (void) malloc_init();
  386.  
  387.     var_label[0] = "x";
  388.     var_label[1] = "y";
  389.     var_label[2] = "t";
  390.     param_label[0] = "omega";
  391.     param_label[1] = "ampl";
  392.     param_label[2] = "damp";
  393.     func_label[0] = "Undefined";
  394.     func_label[1] = "Undefined";
  395.  
  396.     param[0] = 1;
  397.     param[1] = 0;
  398.     param[2] = 0;
  399.  
  400.     var_i[0] = 0;
  401.     var_i[1] = 0;
  402.     var_i[2] = 0;
  403.  
  404.     param_min[0]= -5; param_max[0]= 5;
  405.     param_min[1]= -5; param_max[1]= 5;
  406.     param_min[2]= -5; param_max[2]= 5;
  407.  
  408.     var_min[0]= -5; var_max[0]= 5;
  409.     var_min[1]= -5; var_max[1]= 5;
  410.     var_min[2]= -5; var_max[2]= 5;
  411.  
  412.     f_p = nlmathieu_f;
  413.     func_p = nlmathieu_func;
  414. }
  415. /* nonlinar mathieu equation */
  416. int nlmathieu_f(f,index,x,p,t,dim)
  417. int index,dim;
  418. double f[],x[],p[],t;
  419. {
  420.     f[0] = x[1];
  421.     f[1] = -p[2] * x[1] - (p[0] * p[0] + p[1] * cos(x[2])) * sin(x[0]);
  422.     f[2] = 1;
  423.     
  424. }
  425. int nlmathieu_func(f,x,p,t,dim)
  426. double f[],x[],p[],t;
  427. int dim;
  428. {
  429. }
  430.  
  431. int dpfosc2_init()
  432. {
  433.     title_label = "Diss Period. Forced Oscillator";
  434.  
  435.     mapping_on = 0;
  436.     inverse_on = 1;
  437.     fderiv_on = 0;
  438.     enable_polar = 0;
  439.     enable_period = 1;
  440.     var_dim = 2;
  441.     param_dim = 5;
  442.     func_dim = 2;
  443.  
  444.     (void) malloc_init();
  445.  
  446.     period_len[0] = 1;
  447.  
  448.     var_label[0] = "x";
  449.     var_label[1] = "y";
  450.     param_label[0] = "f-omega";
  451.     param_label[1] = "acampl";
  452.     param_label[2] = "dcampl";
  453.     param_label[3] = "damp";
  454.     param_label[4] = "n-omega";
  455.     func_label[0] = "t";
  456.     func_label[1] = "Undefined";
  457.  
  458.     param[0] = 2;
  459.     param[1] = 1;
  460.     param[2] = 0;
  461.     param[3] = 0;
  462.     param[4] = 0;
  463.  
  464.     var_i[0] = 0;
  465.     var_i[1] = 0.5;
  466.  
  467.     /* stating values of parameter window box */
  468.     param_min[0]= -5; param_max[0]= 5;
  469.     param_min[1]= -5; param_max[1]= 5;
  470.     param_min[2]= -5; param_max[2]= 5;
  471.     param_min[3]= -5; param_max[3]= 5;
  472.     param_min[4]= -5; param_max[4]= 5;
  473.  
  474.     var_min[0]= 0; var_max[0]= 1;
  475.     var_min[1]= -5; var_max[1]= 5;
  476.     func_min[0]= -10; func_max[0]= 10;
  477.  
  478.     f_p = dpfosc2_f;
  479.     func_p = dpfosc2_func;
  480. }
  481. int dpfosc2_f(f,index,x,p,t,dim)
  482. int index,dim;
  483. double f[],x[],p[],t;
  484. {
  485.     f[0] = x[1];
  486.     f[1] = p[2] + p[1]*cos(twopi*p[0]* t) - p[3] * x[1] - p[4] * sin(twopi*x[0]);
  487. }
  488. int dpfosc2_func(f,x,p,t,dim)
  489. double f[],x[],p[],t;
  490. int dim;
  491. {
  492.     f[0] = t;
  493. }
  494. int henonmap_init()
  495. {
  496.     title_label = "Henon Map";
  497.  
  498.     mapping_on = 1;
  499.     inverse_on = 1;
  500.     fderiv_on = 0;
  501.     enable_polar = 0;
  502.     enable_period = 0;
  503.     
  504.     var_dim = 2;
  505.     param_dim = 2;
  506.     func_dim = 0;
  507.  
  508.     (void) malloc_init();
  509.  
  510.     var_label[0] = "x";
  511.     var_label[1] = "y";
  512.     param_label[0] = "a";
  513.     param_label[1] = "b";
  514.  
  515.     param[0] = 1.34;
  516.     param[1] = 0.3;
  517.  
  518.     var_i[0] = 0;
  519.     var_i[1] = 0;
  520.  
  521.     param_min[0]= -5; param_max[0]= 5;
  522.     param_min[1]= -5; param_max[1]= 5;
  523.  
  524.     var_min[0]= -5; var_max[0]= 5;
  525.     var_min[1]= -5; var_max[1]= 5;
  526.  
  527.     f_p = henonmap_f;
  528.     func_p = henonmap_func;
  529. }
  530. /* Henon map */
  531. int henonmap_f(f,index,x,p,t,dim)
  532. int index,dim;
  533. double f[],x[],p[],t;
  534. {
  535.     /* an example of a mapping */
  536.     /* forward mapping */
  537.     if(index ==1) {
  538.         f[0] = 1. + x[1] - p[0] * x[0] * x[0];
  539.         f[1] = p[1] * x[0]; 
  540.     }
  541.     /* backward mapping */
  542.     else if(index ==0) {
  543.         f[0] = 1./p[1] * x[1];
  544.         f[1] = -1. + x[0] + p[0]/p[1]/p[1] * x[1] * x[1]; 
  545.     }
  546.     
  547. }
  548. int henonmap_func(f,x,p,t,dim)
  549. double f[],x[],p[],t;
  550. int dim;
  551. {
  552.     /* an example of no definition */
  553. }
  554. int kotorusmap_init()
  555. {
  556.     title_label = "Kim-Ostlund Torus Map";
  557.     mapping_on = 1;
  558.     inverse_on = 0;
  559.     fderiv_on = 0;
  560.     enable_polar = 0;
  561.     enable_period = 1;
  562.     
  563.     var_dim = 2;
  564.     param_dim = 4;
  565.     func_dim = 2;
  566.  
  567.     (void) malloc_init();
  568.  
  569.     period_len[0] = 1;
  570.     period_len[1] = 1;
  571.     var_label[0] = "x";
  572.     var_label[1] = "y";
  573.     param_label[0] = "wx";
  574.     param_label[1] = "wy";
  575.     param_label[2] = "a";
  576.     param_label[3] = "asymm";
  577.     func_label[0] = "rhox";
  578.     func_label[1] = "rhoy";
  579.  
  580.     param[0] = 0.6;
  581.     param[1] = 0.805;
  582.     param[2] = 0.7;
  583.     param[3] = 1;
  584.     var_i[0] = 0;
  585.     var_i[1] = 0;
  586.  
  587.     param_min[0]= 0; param_max[0]= 1;
  588.     param_min[1]= 0; param_max[1]= 1;
  589.     param_min[2]= 0; param_max[2]= 2;
  590.     param_min[3]= 0; param_max[3]= 2;
  591.     var_min[0]= 0; var_max[0]= 1;
  592.     var_min[1]= 0; var_max[1]= 1;
  593.     func_min[0]= -1; func_max[0]= 1;
  594.     func_min[1]= -1; func_max[1]= 1;
  595.  
  596.     f_p = kotorusmap_f;
  597.     func_p = kotorusmap_func;
  598. }
  599. /* Kim-Ostlund strongly coupled torus map */
  600. int kotorusmap_f(f,index,x,p,t,dim)
  601. int index,dim;
  602. double f[],x[],p[],t;
  603. {
  604.     /* forward map */
  605.     if(index ==1) {
  606.         f[0] = x[0] + p[0] - p[2] * p[3] / twopi * sin(twopi * x[1]);
  607.         f[1] = x[1] + p[1] - p[2] / p[3] / twopi * sin(twopi * x[0]);
  608.     }
  609.     /* inverse map not defined */
  610. }
  611. int kotorusmap_func(f,x,p,t,dim)
  612. double f[],x[],p[],t;
  613. int dim;
  614. {
  615.     int i;
  616.     extern int forward_toggle;
  617.     extern double *v;
  618.     extern int (*f_p)();
  619.     
  620.     (int) f_p(v,forward_toggle,x,p,t,dim);
  621.     f[0] = v[0]-x[0]-p[0];
  622.     f[1] = v[1]-x[1]-p[1];
  623. }
  624. int dissstandmap_init()
  625. {
  626.     title_label = "(Dissipative) Standard Map";
  627.     mapping_on = 1;
  628.     inverse_on = 1;
  629.     fderiv_on = 0;
  630.     enable_polar = 0;
  631.     enable_period = 1;
  632.     var_dim = 2;
  633.     param_dim = 3;
  634.     func_dim = 2;
  635.  
  636.     (void) malloc_init();
  637.  
  638.     period_len[0] = 1;
  639.     period_len[1] = 0;
  640.     
  641.     var_label[0] = "x";
  642.     var_label[1] = "r";
  643.     param_label[0] = "w";
  644.     param_label[1] = "k";
  645.     param_label[2] = "b";
  646.     func_label[0] = "Rhox";
  647.     func_label[1] = "Undefined";
  648.  
  649.     param[0] = 0;
  650.     param[1] = 0.97;
  651.     param[2] = 1;
  652.     var_i[0] = .5;
  653.     var_i[1] = .5;
  654.  
  655.     param_min[0]= 0; param_max[0]= 1;
  656.     param_min[1]= 0; param_max[1]= 1;
  657.     param_min[2]= 0; param_max[2]= 1;
  658.     var_min[0]= 0; var_max[0]= 1;
  659.     var_min[1]= 0; var_max[1]= 1;
  660.     func_min[0]= 0; func_max[0]= 1;
  661.  
  662.     f_p = dissstandmap_f;
  663.     func_p = dissstandmap_func;
  664. }
  665. /* Diss. Standard Map */
  666. int dissstandmap_f(f,index,x,p,t,dim)
  667. int index,dim;
  668. double f[],x[],p[],t;
  669. {
  670.     double v0sq,v1sq,v2sq,v3sq;
  671.  
  672.     /* forward map */
  673.     if(index ==1) {
  674.         f[1] = p[2] * x[1] - p[1] / twopi * sin(twopi * x[0]);
  675.         f[0] = x[0] + p[0] + f[1];
  676.     }
  677.     /* backward map */
  678.     else if(index ==0) {
  679.         f[0] = x[0] - p[0] -x[1];
  680.         f[1] = (x[1] + p[1]/twopi * sin(twopi * f[0]))/p[2];
  681.     }
  682.     
  683. }
  684. int dissstandmap_func(f,x,p,t,dim)
  685. double f[],x[],p[],t;
  686. int dim;
  687. {
  688.     int i;
  689.     extern int forward_toggle;
  690.     extern double *v;
  691.     extern int (*f_p)();
  692.     
  693.     (int) f_p(v,forward_toggle,x,p,t,dim);
  694.     f[0] = v[0]-x[0];
  695. }
  696. int siegelmap_init()
  697. {
  698.     title_label = "Siegel Map";
  699.  
  700.     mapping_on = 1;
  701.     inverse_on = 0;
  702.     fderiv_on = 0;
  703.     enable_polar = 0;
  704.     enable_period = 0;
  705.     
  706.     var_dim = 2;
  707.     param_dim = 2;
  708.     func_dim = 2;
  709.  
  710.     (void) malloc_init();
  711.  
  712.     var_label[0] = "x";
  713.     var_label[1] = "y";
  714.     param_label[0] = "rho";
  715.     param_label[1] = "p";
  716.     func_label[0] = "|Rho|^2";
  717.     func_label[1] = "Undefined";
  718.  
  719.     param[0] = 0.618;
  720.     param[1] = 1.;
  721.     var_i[0] = 1.;
  722.     var_i[1] = 0.;
  723.  
  724.     param_min[0]= 0; param_max[0]= 1;
  725.     param_min[1]= 0; param_max[1]= 10;
  726.     var_min[0]= -2; var_max[0]= 2;
  727.     var_min[1]= -2; var_max[1]= 2;
  728.     func_min[0]= 0; func_max[0]= 4;
  729.  
  730.     f_p = siegelmap_f;
  731.     func_p = siegelmap_func;
  732. }
  733. /* Siegel map */
  734. int siegelmap_f(f,index,x,p,t,dim)
  735. int index,dim;
  736. double f[],x[],p[],t;
  737. {
  738.     double xp,yp,xt,yt,theta,rr,alpha,sigma,cx,cy,cossigma,sinsigma;
  739.     if(index ==1) {
  740.         xt = x[0];
  741.         yt= x[1];
  742.         alpha = p[1];
  743.         sigma = p[0];
  744.         sinsigma = sin(twopi*sigma);
  745.         cossigma = cos(twopi*sigma);
  746.         xp = 1-xt;
  747.         yp = -yt;
  748.         if(xp==0){
  749.             if(yp>0)
  750.                 theta = twopi/4;
  751.             else
  752.                 theta = -twopi/4;
  753.         }
  754.         else {
  755.             theta = atan(yp/xp);
  756.             if(xp>0 && yp>=0)
  757.                 theta = theta; 
  758.             else if(xp>0 && yp<0)
  759.                 theta = theta;
  760.             else if(xp<0 && yp>0)
  761.                 theta = twopi/2 + theta;
  762.             else if(xp<0 && yp<0)
  763.                 theta = -(twopi/2 - theta);
  764.         }
  765.         rr = xp*xp + yp*yp;
  766.         rr = exp((alpha+1)/2 * log(rr));
  767.         cx = 1 - rr * cos((alpha+1) * theta);
  768.         cy = -rr * sin((alpha+1) * theta);
  769.         cx = cx / (alpha+1);
  770.         cy = cy / (alpha+1);
  771.         f[0] = cossigma * cx - sinsigma * cy; 
  772.         f[1] = sinsigma * cx + cossigma * cy; 
  773.     }
  774.     /* inverse map not defined */
  775. }
  776. int siegelmap_func(f,x,p,t,dim)
  777. double f[],x[],p[],t;
  778. int dim;
  779. {
  780.     int i;
  781.     extern int forward_toggle;
  782.     extern double *v;
  783.     extern int (*f_p)();
  784.     
  785.     (int) f_p(v,forward_toggle,x,p,t,dim);
  786.     f[0] = 0;
  787.     for(i=0;i<dim;i++)
  788.         f[0] += (v[i]-x[i])*(v[i]-x[i]);
  789.     
  790. }
  791. /*
  792. This is a martyd3 subroutine. "martyd3" should be
  793. substituted with the proper string for a dynamical system
  794. to be installed. Define all initialization parameters
  795. here. Parameters are assigned to the default values before this program is
  796. called.
  797. */
  798. int martyd3_init()
  799. {
  800.     title_label = "Marty's D3 Map";
  801.  
  802.     mapping_on = 1;
  803.     inverse_on = 0;
  804.     fderiv_on = 0;
  805.     enable_polar = 0;
  806.     enable_period = 0;
  807.  
  808.     var_dim = 2;
  809.     param_dim = 4;
  810.     func_dim = 2;
  811.  
  812.     (void) malloc_init();
  813.  
  814.     var_label[0] = "x";
  815.     var_label[1] = "y";
  816.     param_label[0] = "alpha";
  817.     param_label[1] = "lambda";
  818.     param_label[2] = "beta";
  819.     param_label[3] = "gamma";
  820.     func_label[0] = "Modulus";
  821.     func_label[1] = "Re(Z^3)";
  822.  
  823.     param[0] = -1;
  824.     param[1] = 1.52;
  825.     param[2] = .1;
  826.     param[3] = -.8;
  827.     var_i[0] = 0.01;
  828.     var_i[1] = 0.057;
  829.  
  830.     param_min[0]= -5; param_max[0]= 5;
  831.     param_min[1]= -5; param_max[1]= 5;
  832.     param_min[2]= -5; param_max[2]= 5;
  833.     param_min[3]= -5; param_max[3]= 5;
  834.     var_min[0]= -1.5; var_max[0]= 1.5;
  835.     var_min[1]= -1.5; var_max[1]= 1.5;
  836.  
  837.     f_p = martyd3_f;
  838.     func_p = martyd3_func;
  839. }
  840. /*
  841. second user dynamical system 
  842. */
  843.     
  844. int martyd3_f(f,index,x,p,t,dim)
  845. int index,dim;
  846. double f[],x[],p[],t;
  847. {
  848.     double z2r,z2i,z3r,z3i,iv;
  849.         cmul(&z2r,&z2i,x[0],x[1],x[0],x[1]);
  850.     cmul(&z3r,&z3i,z2r,z2i,x[0],x[1]);
  851.     iv = p[0] * (x[0]*x[0] + x[1]*x[1]) + p[1] + p[2] * z3r; 
  852.     if(index ==1) {
  853.         f[0] = iv * x[0] + p[3] * z2r;
  854.         f[1] = iv * x[1] - p[3] * z2i;
  855.     }
  856. }
  857.  
  858. cmul(x,y,x1,y1,x2,y2)
  859. double *x,*y,x1,y1,x2,y2;
  860. {
  861.     *x = x1 * x2 - y1 * y2; 
  862.     *y = x1 * y2 + x2 * y1; 
  863. }
  864. /* cmul() is a subroutine local for this model */
  865. int martyd3_func(f,x,p,t,dim)
  866. double f[],x[],p[],t;
  867. int dim;
  868. {
  869.     double z2r,z2i,z3r,z3i,iv;
  870.         cmul(&z2r,&z2i,x[0],x[1],x[0],x[1]);
  871.     cmul(&z3r,&z3i,z2r,z2i,x[0],x[1]);
  872.     f[0] = x[0]*x[0] + x[1]*x[1];
  873.     f[1] = z3r;
  874. }
  875.  
  876. int henonheiles_init()
  877. {
  878.     title_label = "Henon-Heiles";
  879.  
  880.     mapping_on = 0;
  881.     inverse_on = 1;
  882.     fderiv_on = 0;
  883.     enable_polar = 0;
  884.     enable_period = 0;
  885.     var_dim = 4;
  886.     param_dim = 1;
  887.     func_dim = 3;
  888.  
  889.     (void) malloc_init();
  890.  
  891.     var_label[0] = "x";
  892.     var_label[1] = "px";
  893.     var_label[2] = "y";
  894.     var_label[3] = "py";
  895.     param_label[0] = "epsilon";
  896.     func_label[0] = "Energy";
  897.     func_label[1] = "Ang Mon";
  898.     func_label[2] = "t";
  899.  
  900.     param[0] = 1;
  901.  
  902.     var_i[0] = 0.1;
  903.     var_i[1] = 0.1;
  904.     var_i[2] = 0.1;
  905.     var_i[3] = 0.1;
  906.  
  907.     /* stating values of parameter window box */
  908.     param_min[0]= 0; param_max[0]= 5;
  909.  
  910.     var_min[0]= -1; var_max[0]= 1;
  911.     var_min[1]= -1; var_max[1]= 1;
  912.     var_min[2]= -1; var_max[3]= 1;
  913.     var_min[3]= -1; var_max[2]= 1;
  914.     func_min[0]= -1; func_max[0]= 1;
  915.     func_min[1]= -1; func_max[1]= 1;
  916.     func_min[2]= 0; func_max[2]= 10000;
  917.  
  918.     f_p = henonheiles_f;
  919.     func_p = henonheiles_func;
  920. }
  921.  
  922. int henonheiles_f(f,index,x,p,t,dim)
  923. int index,dim;
  924. double f[],x[],p[],t;
  925. {
  926.     if(index!=2){
  927.         f[0] = x[1];
  928.         f[2] = x[3];
  929.     }
  930.     if(index!=1){
  931.         f[1] = - x[0] - p[0]*2.*x[0]*x[2];
  932.         f[3] = - x[2] - p[0]*(x[0]*x[0]-x[2]*x[2]);
  933.     }
  934. }
  935.  
  936. int henonheiles_func(f,x,p,t,dim)
  937. double f[],x[],p[],t;
  938. int dim;
  939. {
  940.     double x0sq,x2sq;
  941.     x0sq = x[0]*x[0];
  942.     x2sq = x[2]*x[2];
  943.     f[0] = 0.5 * (x[1]*x[1] + x[3]*x[3] + x0sq + x2sq) + p[0]*(x0sq*x[2] - x2sq*x[2]/3.);
  944.     f[1] = x[0]*x[3] - x[1]*x[2];
  945.     f[2] = t;
  946. }
  947. int vanderpol_init()
  948. {
  949.     title_label = "Van der Pol's Eq";
  950.  
  951.     mapping_on = 0;
  952.     inverse_on = 1;
  953.     fderiv_on = 0;
  954.     enable_polar = 0;
  955.     enable_period = 0;
  956.  
  957.     var_dim = 2;
  958.     param_dim = 3;
  959.     func_dim = 2;
  960.  
  961.     (void) malloc_init();
  962.  
  963.     var_label[0] = "x";
  964.     var_label[1] = "y";
  965.     param_label[0] = "alpha";
  966.     param_label[1] = "beta";
  967.     param_label[2] = "omega";
  968.     func_label[0] = "t";
  969.     func_label[1] = "Undefined";
  970.  
  971.     param[0] = 1;
  972.     param[1] = 0;
  973.     param[2] = 1;
  974.  
  975.     var_i[0] = 0;
  976.     var_i[1] = 0;
  977.  
  978.     param_min[0]= -5; param_max[0]= 5;
  979.     param_min[1]= -5; param_max[1]= 5;
  980.     param_min[2]= 0; param_max[2]= 5;
  981.  
  982.     var_min[0]= -5; var_max[0]= 5;
  983.     var_min[1]= -5; var_max[1]= 5;
  984.  
  985.     f_p = vanderpol_f;
  986.     func_p = vanderpol_func;
  987. }
  988. /* Van der Pol's equation */
  989. int vanderpol_f(f,index,x,p,t,dim)
  990. int index,dim;
  991. double f[],x[],p[],t;
  992. {
  993.     double cos();
  994.     f[0] = x[1] - p[0]*(x[0]*x[0]*x[0]/3. - x[0]);
  995.     f[1] = -x[0] + p[1]*cos(p[2]*t);
  996. }
  997. int vanderpol_func(f,x,p,t,dim)
  998. double f[],x[],p[],t;
  999. int dim;
  1000. {
  1001.     f[0] = t;
  1002.     f[1] = 0;
  1003. }
  1004. int duffing_init()
  1005. {
  1006.     title_label = "Duffing's Eq";
  1007.  
  1008.     mapping_on = 0;
  1009.     inverse_on = 1;
  1010.     fderiv_on = 0;
  1011.     enable_polar = 0;
  1012.     enable_period = 0;
  1013.  
  1014.     var_dim = 2;
  1015.     param_dim = 4;
  1016.     func_dim = 2;
  1017.  
  1018.     (void) malloc_init();
  1019.  
  1020.     var_label[0] = "x";
  1021.     var_label[1] = "p";
  1022.     param_label[0] = "delta";
  1023.     param_label[1] = "beta";
  1024.     param_label[2] = "gamma";
  1025.     param_label[3] = "omega";
  1026.     func_label[0] = "t";
  1027.     func_label[1] = "Undefined";
  1028.  
  1029.     param[0] = 0.25;
  1030.     param[1] = 1;
  1031.     param[2] = 0.3;
  1032.     param[3] = 1;
  1033.  
  1034.     var_i[0] = 0;
  1035.     var_i[1] = 0;
  1036.  
  1037.     param_min[0]= -5; param_max[0]= 5;
  1038.     param_min[1]= -5; param_max[1]= 5;
  1039.     param_min[2]= -5; param_max[2]= 5;
  1040.     param_min[3]= 0; param_max[3]= 5;
  1041.  
  1042.     var_min[0]= -5; var_max[0]= 5;
  1043.     var_min[1]= -5; var_max[1]= 5;
  1044.  
  1045.     f_p = duffing_f;
  1046.     func_p = duffing_func;
  1047. }
  1048. /* Duffing's equation */
  1049. int duffing_f(f,index,x,p,t,dim)
  1050. int index,dim;
  1051. double f[],x[],p[],t;
  1052. {
  1053.     double cos();
  1054.     f[0] = x[1];
  1055.     f[1] = p[1]*x[0] - x[0]*x[0]*x[0] - p[0]*x[1] + p[2]*cos(p[3]*t);
  1056. }
  1057. int duffing_func(f,x,p,t,dim)
  1058. double f[],x[],p[],t;
  1059. int dim;
  1060. {
  1061.     f[0] = t;
  1062.     f[1] = 0;
  1063. }
  1064. int simpletorusmap_init()
  1065. {
  1066.     title_label = "Simple Torus Map";
  1067.     mapping_on = 1;
  1068.     inverse_on = 0;
  1069.     fderiv_on = 0;
  1070.     enable_polar = 0;
  1071.     enable_period = 1;
  1072.     
  1073.     var_dim = 2;
  1074.     param_dim = 4;
  1075.     func_dim = 2;
  1076.  
  1077.     (void) malloc_init();
  1078.  
  1079.     period_len[0] = 1;
  1080.     period_len[1] = 1;
  1081.     var_label[0] = "x";
  1082.     var_label[1] = "y";
  1083.     param_label[0] = "wx";
  1084.     param_label[1] = "wy";
  1085.     param_label[2] = "a";
  1086.     param_label[3] = "e";
  1087.     func_label[0] = "rhox";
  1088.     func_label[1] = "rhoy";
  1089.  
  1090.     param[0] = 0.1;
  1091.     param[1] = 1.;
  1092.     param[2] = 0.2;
  1093.     param[3] = 0.1;
  1094.     var_i[0] = 0;
  1095.     var_i[1] = 0;
  1096.  
  1097.     param_min[0]= -1.5; param_max[0]= 1.5;
  1098.     param_min[1]= -1.5; param_max[1]= 1.5;
  1099.     param_min[2]= 0; param_max[2]= 2;
  1100.     param_min[3]= 0; param_max[3]= 2;
  1101.     var_min[0]= 0; var_max[0]= 1;
  1102.     var_min[1]= 0; var_max[1]= 1;
  1103.     func_min[0]= -1; func_max[0]= 1;
  1104.     func_min[1]= -1; func_max[1]= 1;
  1105.  
  1106.     f_p = simpletorusmap_f;
  1107.     func_p = simpletorusmap_func;
  1108. }
  1109. /* Kim-Ostlund strongly coupled torus map */
  1110. int simpletorusmap_f(f,index,x,p,t,dim)
  1111. int index,dim;
  1112. double f[],x[],p[],t;
  1113. {
  1114.     /* forward map */
  1115.     if(index ==1) {
  1116.         f[0] = x[0] + p[3]*(p[0] + cos(2*pi*x[0]) + p[2]*cos(2*pi*x[1]));
  1117.         f[1] = x[1] + p[3]*(p[1] + sin(2*pi*x[0]) + p[2]*sin(2*pi*x[1]));
  1118.     }
  1119.     /* inverse map not defined */
  1120. }
  1121. int simpletorusmap_func(f,x,p,t,dim)
  1122. double f[],x[],p[],t;
  1123. int dim;
  1124. {
  1125.     int i;
  1126.     extern int forward_toggle;
  1127.     extern double *v;
  1128.     extern int (*f_p)();
  1129.     
  1130.     (int) f_p(v,forward_toggle,x,p,t,dim);
  1131.     f[0] = v[0]-x[0];
  1132.     f[1] = v[1]-x[1];
  1133. }
  1134.