home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / modellib / generic_def.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-16  |  11.7 KB  |  431 lines
  1. /*
  2. Initialize all parameters for a given dynamical system to be installed
  3. Parameters are assigned to the default values before this program is called.
  4. -----------------------------------------------------------------------
  5. This is a GENERIC subroutine. The string "generic" should be
  6. GLOBALLY substituted with the proper name string for a given dynamical system
  7. to be installed.
  8. */
  9.  
  10. /*
  11. Example 1: 4-d vector field with polar coords and no periodic variable
  12. */
  13. int generic_init()
  14. {
  15.     /* title label */
  16.     title_label = "Undefined GENERIC case";
  17.  
  18.     /* mapping toggle: 1=map, 0=vector field */
  19.     mapping_on = 0;
  20.     /* inverse toggle: vector field: always 1,
  21.        maps: 0=explicit inverse is defined, 1=otherwise  */
  22.     inverse_on = 1;
  23.     /* jacobian toggle: 1=Jacobian is explicitly given, 0=otherwise */
  24.     fderiv_on = 0;
  25.     /* polar toggle: 1=enable polar coordinate feature, 0=otherwise */
  26.     enable_polar = 1;
  27.     /* period toggle: 1=enable periodicity of phase space, 0=otherwise */
  28.     enable_period = 0;
  29.  
  30.     /* phase space dimension */
  31.     var_dim = 4;
  32.     /* parameter space dimension */
  33.     param_dim = 2;
  34.     /* function space dimension */
  35.     func_dim = 2;
  36.  
  37.     (void) malloc_init();
  38.  
  39.     /*periodicity of phase space variables (DIM=var_dim)*/
  40.     /* do not have to be defined if enable_period = 0 */
  41.     period_len[0] =0;
  42.     period_len[1] =0;
  43.     period_len[2] =0;
  44.     period_len[3] =0;
  45.  
  46.     /* primary phase space variable label (DIM=var_dim)*/
  47.     var_label[0] = "x";
  48.     var_label[1] = "y";
  49.     var_label[2] = "z";
  50.     var_label[3] = "w";
  51.     /* secondary phase space variable label (usually polar coord variable)
  52.     (DIM=var_dim) */
  53.     /* do not have to be defined if enable_polar = 0 */
  54.     var_polar_label[0] = "r";
  55.     var_polar_label[1] = "rp";
  56.     var_polar_label[2] = "theta";
  57.     var_polar_label[3] = "thetap";
  58.     /* parameter variable label (DIM=param_dim)*/
  59.     param_label[0] = "mu";
  60.     param_label[1] = "delta";
  61.     /* function variable label (DIM=func_dim)*/
  62.     func_label[0] = "Energy";
  63.     func_label[1] = "AngMom";
  64.  
  65.     /* starting parameter values (DIM=param_dim)*/
  66.     param[0] = 2;
  67.     param[1] = 1;
  68.  
  69.     /* starting primary phase space variable values (DIM=param_dim)*/
  70.     var_i[0] = 0;
  71.     var_i[1] = 0;
  72.     var_i[2] = 0;
  73.     var_i[3] = 0;
  74.     /* starting seconsary phase space variable values (DIM=param_dim)*/
  75.     /* do not have to be defined if enable_polar = 0 */
  76.     var_polar_i[0] = 0;
  77.     var_polar_i[1] = 0;
  78.     var_polar_i[2] = 0;
  79.     var_polar_i[3] = 0;
  80.  
  81.     /* starting bounds of parameter window box */
  82.     param_min[0]= -5; param_max[0]= 5;
  83.     param_min[1]= -5; param_max[1]= 5;
  84.  
  85.     /* starting bounds of primary phase space window box */
  86.     var_min[0]= -5; var_max[0]= 5;
  87.     var_min[1]= -5; var_max[1]= 5;
  88.     var_min[2]= -5; var_max[2]= 5;
  89.     var_min[3]= -5; var_max[3]= 5;
  90.  
  91.     /* starting bounds of secondary phase space window box */
  92.     /* do not have to be defined if enable_polar = 0 */
  93.     var_polar_min[0]= -5; var_polar_max[0]= 5;
  94.     var_polar_min[1]= -5; var_polar_max[1]= 5;
  95.     var_polar_min[2]= -pi; var_polar_max[2]= pi;
  96.     var_polar_min[3]= -5; var_polar_max[3]= 5;
  97.  
  98.     /* dynamical system and function pointer assignments */
  99.     f_p = generic_f;
  100.     func_p = generic_func;
  101. }
  102.  
  103. /*
  104. Example 2: 2-d torus map with periodic variables and no polar coords  
  105. */
  106. int userds0_init()
  107. {
  108.     extern int (*f_p)(),userds0_f(),(*func_p)(),userds0_func();
  109.     extern int var_dim,func_dim,param_dim;
  110.     extern int mapping_on,inverse_on,fderiv_on,enable_polar,enable_period;
  111.     extern double *param,*param_min,*param_max,*func,*func_min,*func_max;
  112.     extern double *var_i,*var_polar_i,*var_min,*var_max,*var_polar_min,*var_polar_max;
  113.     extern double *period_len,pi;
  114.     extern char *title_label,**var_label,**var_polar_label,**param_label,**func_label;
  115.  
  116.     /* title label */
  117.     title_label = "Undefined GENERIC case";
  118.  
  119.     /* mapping toggle: 1=map 0=vector field */
  120.     mapping_on = 1;
  121.     /* inverse toggle: vector field: always 1,
  122.        maps: 0=explicit inverse is defined, 1=otherwise  */
  123.     inverse_on = 0;
  124.     /* jacobian toggle: 1=Jacobian is explicitly given, 0=otherwise */
  125.     fderiv_on = 0;
  126.     /* polar toggle: 1=enable polar coordinate feature, 0=otherwise */
  127.     enable_polar = 0;
  128.     /* period toggle: 1=enable periodicity of phase space, 0=otherwise */
  129.     enable_period = 1;
  130.  
  131.     /* phase space dimension */
  132.     var_dim = 2;
  133.     /* parameter space dimension */
  134.     param_dim = 2;
  135.  
  136.     /*periodicity of phase space variables (DIM=var_dim)*/
  137.     /* do not have to be defined if enable_period = 0 */
  138.     period_len[0] =1;
  139.     period_len[1] =1;
  140.  
  141.     /* primary phase space variable label (DIM=var_dim)*/
  142.     var_label[0] = "x";
  143.     var_label[1] = "y";
  144.     /* seconsary phase space variable label (DIM=param_dim)*/
  145.     /* do not have to be defined if enable_polar = 0 */
  146.     
  147.     /* parameter variable label (DIM=param_dim)*/
  148.     param_label[0] = "a";
  149.     param_label[1] = "b";
  150.     /* function variable label (DIM=func_dim)*/
  151.     func_label[0] = "F1";
  152.     func_label[1] = "F2";
  153.  
  154.     /* starting parameter values (DIM=param_dim)*/
  155.     param[0] = 2;
  156.     param[1] = 1;
  157.  
  158.     /* starting primary phase space variable values (DIM=param_dim)*/
  159.     var_i[0] = 0;
  160.     var_i[1] = 0;
  161.     /* starting seconsary phase space variable values (DIM=param_dim)*/
  162.     /* do not have to be defined if enable_polar = 0 */
  163.  
  164.     /* starting bounds of parameter window box */
  165.     param_min[0]= -5; param_max[0]= 5;
  166.     param_min[1]= -5; param_max[1]= 5;
  167.  
  168.     /* starting bounds of primary phase space window box */
  169.     var_min[0]= 0; var_max[0]= 1;
  170.     var_min[1]= 0; var_max[1]= 1;
  171.  
  172.     /* starting bounds of secondary phase space window box */
  173.     /* do not have to be defined if enable_polar = 0 */
  174.  
  175.     /* dynamical system and function pointer assignments */
  176.     f_p = userds0_f;
  177.     func_p = userds0_func;
  178. }
  179. /*
  180. -----------------------------------------------------------------
  181. user definition of an ordinary differential equation or a mapping
  182. -----------------------------------------------------------------
  183. FILES NEED TO CHANGED: generic_func.c
  184.     (optional) generic_choose_algorithm.c
  185. -----------------------------------------------------------------
  186. Input:    x[n]: coords
  187.     p[n]: parameters
  188.     t: time (used for ode only)
  189.     dim: dimension of coords
  190.     index:    (ode) 0: all equations
  191.         (Hamiltonian ode only) 1: momentum part, 2: coord part
  192.         (map) 1: forward, 0: backward
  193. Output:    f[n]: (ode) values of vector fields, (map) an image
  194.     (values of the right hand side of dynamical system equation)
  195. -----------------------------------------------------------------
  196. Useful system parameters which have been already declared and can be used here:
  197. -----------------------------------------------------------------
  198.     (double) pi:    Pi
  199.     (double) twopi:    2 * Pi
  200. -----------------------------------------------------------------
  201. Useful system parameters which can be externally declared and used here:
  202. -----------------------------------------------------------------
  203.     (int) var_dim: dimension of the phase space
  204.     (int) param_dim: dimension of the system parameters
  205.     (int) model_dim: dimension of a list of all models installed
  206.     (int) model: currently selected model
  207.     (int) cur_color: currently selected color (sometimes overriden)
  208.     (int) symbol_type: currently selected symbol (sometimes overriden)
  209.     (int) symbol_size: size of an symbol (sometimes overriden)
  210. -----------------------------------------------------------------
  211. */
  212.     
  213. /*
  214. Example 1: mapping without an explicit inverse
  215. */
  216. int generic_f(f,index,x,p,t,dim)
  217. int index,dim;
  218. double f[],x[],p[],t;
  219. {
  220.     f[0] = p[0] * ( x[1] - x[0] );
  221.     f[1] = p[1] * x[0] - x[1] - x[0] * x[2];
  222.     f[2] = - p[2] * x[2] + x[0] * x[1];
  223. }
  224.  
  225. /*
  226. Example 2: ODE (non-Hamiltonian)
  227. Note that there is an explicit dependence on time t
  228. and sin and cos are declared as double
  229. */
  230. int generic_f(f,index,x,p,t,dim)
  231. int index,dim;
  232. double f[],x[],p[],t;
  233. {
  234.     f[0] = x[1];
  235.     f[1] = -p[2] * x[1] - (p[0] * p[0] + p[1] * cos(t) * sin(x[0]);
  236. }
  237.  
  238. /*
  239. Example 3: mapping with an explicit inverse
  240. */
  241. int generic_f(f,index,x,p,t,dim)
  242. int index,dim;
  243. double f[],x[],p[],t;
  244. {
  245.     if(index ==1) {
  246.         /* forward mapping */
  247.         f[0] = 1. + x[1] - p[0] * x[0] * x[0];
  248.         f[1] = p[1] * x[0]; 
  249.     }
  250.     else if(index == 0){
  251.         /* backward mapping */
  252.         f[0] = 1. / p[1] * x[1];
  253.         f[1] = -1. + x[0] + p[0] / p[1] /p[1] * x[1] * x[1]; 
  254.     }
  255. }
  256.  
  257. /*
  258. Example 4: Hamiltonian vector field (when symplectic integration is not used)
  259. Note: Often introduction of new variables reduces the number of
  260. multiplcations, and therefore computation time
  261. */
  262. int generic_f(f,index,x,p,t,dim)
  263. int index,dim;
  264. double f[],x[],p[],t;
  265. {
  266.     double v0sq,v2sq;
  267.  
  268.     v0sq = x[0] * x[0];
  269.     v2sq = x[2] * x[2];
  270.     f[0] = x[1];
  271.     f[2] = x[3];
  272.     f[1] = x[0] * (p[0] - (v0sq + v2sq)) + p[1] * x[0] * v2sq;
  273.     f[3] = x[2] * (p[0] - (v0sq + v2sq)) + p[1] * x[2] * v0sq; 
  274. }
  275.  
  276. /*
  277. Example 5: Hamiltonian vector field (when symplectic integration is used)
  278. it should be divided into two sets of equations,
  279. momentum part: dx/dt = dH/dp, coordinate part: dp/dt = -dH/dx
  280. */
  281. int generic_f(f,index,x,p,t,dim)
  282. int index,dim;
  283. double f[],x[],p[],t;
  284. {
  285.     double v0sq,v2sq;
  286.  
  287.     if(index !=2) {
  288.         /* momenturm part dH/dp */
  289.         f[0] = x[1];
  290.         f[2] = x[3];
  291.     }
  292.     if(index !=1){
  293.         /* coordinate part -dH/dx */
  294.         v0sq = x[0] * x[0];
  295.         v2sq = x[2] * x[2];
  296.         f[1] = x[0] * (p[0] - (v0sq + v2sq)) + p[1] * x[0] * v2sq;
  297.         f[3] = x[2] * (p[0] - (v0sq + v2sq)) + p[1] * x[2] * v0sq; 
  298.     }
  299. }
  300.  
  301. /*
  302. Example 6: Mapping without an explicit inverse (use of pi,twopi)
  303. */
  304. int generic_f(f,index,x,p,t,dim)
  305. int index,dim;
  306. double f[],x[],p[],t;
  307. {
  308.     /* forward map */
  309.     f[0] = x[0] + p[0] - p[2] * p[3] / twopi * sin(twopi * x[1]);
  310.     f[1] = x[1] + p[1] - p[2] / p[3] / twopi * sin(twopi * x[0]);
  311. }
  312. /*
  313. Definitions of functions
  314. Use:    (1)fuctions can be regarded as parts of generalized coordinates,
  315.     allowing them to be manipulated as other coordinates.
  316.     For example, one can plot energy versus angular momentum.
  317.     (2)values of functions can be monitored and shown on panels 
  318.     (3)functions can be used to define a general surface for a Poincare
  319.     section.
  320.     (4)functions can be used to define coordinate transformations
  321.     such as polar coordinates (this feature is duplicated by a set
  322.     of predefined coordinate transformation subroutines)
  323.     (5)functions can be colorcoded according to their values 
  324. Note:    functions should be defined in terms of primary coordinates,
  325.     the coordinate system with which you define the dynamical system
  326. */
  327.  
  328. /*
  329. Example 1: no definition
  330. */
  331. int generic_func(f,x,p,t,dim)
  332. double f[],x[],t,p[];
  333. int dim;
  334. {
  335. }
  336.  
  337. /*
  338. Example 2: Energy and angular momentum
  339. */
  340. int generic_func(f,x,p,t,dim)
  341. double f[],x[],p[],t;
  342. int dim;
  343. {
  344.         double v0sq,v1sq,v2sq,v3sq;
  345.  
  346.     coordinates */
  347.     v0sq = x[0] * x[0];
  348.     v2sq = x[2] * x[2];
  349.     f[0] = 0.5 * (( x[1] * x[1] + x[3] * x[3]) - p[0]*(v0sq+v2sq) + 0.5 *
  350.     (v0sq+v2sq)*(v0sq+v2sq) - p[1] * v0sq* v2sq);
  351.     f[1] =  x[1] * x[2] - x[0] * x[3];
  352. }
  353.  
  354. /*
  355. Example 3: in polar coordinates
  356. */
  357. int func_1(f,x,p,t,dim)
  358. double f[],x[],t,p[];
  359. int dim;
  360. {
  361.     double v0sq,v1sq,v2sq,v3sq,sin(),cos(),vsin,vcos;
  362.     extern double *v;
  363.  
  364.     /* polar coordinate to euclidean coordinate transformation */
  365.     euclid_to_polar(v,x);
  366.     
  367.     v0sq = v[0] * v[0];
  368.     v1sq = v[1] * v[1];
  369.     v2sq = v[2] * v[2];
  370.     v3sq = v[3] * v[3];
  371.     vsin = sin(v[2]);
  372.     vcos = cos(v[2]);
  373.  
  374.     /* functions are defined in terms of polar coordinates */
  375.     f[0] = 0.5 * ((v1sq + v0sq * v3sq) - p[0]*v0sq + 0.5 * v0sq * v0sq - p[1] * v0sq * v0sq * vsin * vsin * vcos * vcos);
  376.     f[1] =  v0sq * v[3];
  377. }
  378.  
  379. /*
  380. Example 4: more sophistigated use
  381. */
  382. int generic_func(f,x,p,t,dim)
  383. double f[],x[],t,p[];
  384. int dim;
  385. {
  386.     int i;
  387.     extern int forward_toggle;
  388.     extern double *v;
  389.     extern int (*f_p)();
  390.     
  391.     (int) f_p(v,forward_toggle,x,p,t,dim);
  392.     f[0] = 0;
  393.     for(i=0;i<dim;i++)
  394.         f[0] += (v[i]-x[i])*(v[i]-x[i]);
  395. }
  396.  
  397. /*
  398. Example 5: using time as a function (useful for a conventional time series plot)
  399. */
  400. int generic_func(f,x,p,t,dim)
  401. double f[],x[],t,p[];
  402. int dim;
  403. {
  404.     f[0] = t;
  405. }
  406.  
  407. /*
  408. Example 6: motinor a modulus
  409. */
  410. int generic_func(f,x,p,t,dim)
  411. double f[],x[],t,p[];
  412. int dim;
  413. {
  414.     f[0] = x[0]*x[0] + x[1]*x[1];
  415. }
  416.  
  417. /*
  418. Example 7: partial polar coordinate transformation 
  419. NOTE: By default, some coord transformations are supported by the system.
  420. */
  421. int generic_func(f,x,p,t,dim)
  422. double f[],x[],t,p[];
  423. int dim;
  424. {
  425.     f[0] = x[0]*x[0] + x[1]*x[1];
  426.     if(x[0] != 0)
  427.         f[1] = atan(x[1]/x[0]);
  428.     else
  429.         f[1] = twopi/4.;
  430. }
  431.