PaGMO  1.1.5
taylor_fixedthrust.cpp
1 #include"taylor_fixedthrust.h"
2 
3 
4 int taylor_fixedthrust(double y[6], double t0, double tf,double thrust[3],double abstol, double reltol){
5  int i, j, order=20, itmp=0, direction = 1, nsteps = -1;
6  double ftmp, dstep, tolerance, rtolerance, log10tolerance, log10rtolerance;
7  MY_FLOAT startT, stopT, nextT;
8  MY_FLOAT myFloatZero;
9  MY_FLOAT xx[10], yy[10], zz[10];
10  MY_FLOAT **jet;
11  InitMyFloat(myFloatZero);
12  InitMyFloat(startT);
13  InitMyFloat(stopT);
14  InitMyFloat(nextT);
15  for(i=0; i<10; i++){InitMyFloat(xx[i]); InitMyFloat(yy[i]);InitMyFloat(zz[i]);}
16 
17  MakeMyFloatA(myFloatZero, 0);
18 
19  /* assign initials */
20  MakeMyFloatA(xx[0], (double)y[0] );
21  MakeMyFloatA(xx[1], (double)y[1] );
22  MakeMyFloatA(xx[2], (double)y[2] );
23  MakeMyFloatA(xx[3], (double)y[3] );
24  MakeMyFloatA(xx[4], (double)y[4] );
25  MakeMyFloatA(xx[5], (double)y[5] );
26  MakeMyFloatA(xx[6], (double)thrust[0] );
27  MakeMyFloatA(xx[7], (double)thrust[1] );
28  MakeMyFloatA(xx[8], (double)thrust[2] );
29  MakeMyFloatA(startT, t0);
30  MakeMyFloatA(stopT, tf);
31  dstep=0.001; /* only nedeed when step_ctrl_method==0 (see manual) */
32  MakeMyFloatA(nextT, (double)dstep);
33 
34  log10tolerance = log10(abstol);
35  log10rtolerance = log10(reltol);
36 
37  /* the main loop */
38 
39  if(dstep < (double)0.0) { direction = -1;}
40  do {
41 
42  itmp = taylor_step_kepler_info( &startT, xx, direction, 1, log10tolerance, log10rtolerance, &stopT, &nextT, &order);
43 
44  } while(itmp == 0); /* while */
45  y[0]=xx[0];
46  y[1]=xx[1];
47  y[2]=xx[2];
48  y[3]=xx[3];
49  y[4]=xx[4];
50  y[5]=xx[5];
51 
52  // TODO: check if this return value has meaning or purpose.
53  return 0;
54 }
55 
56 
57 
58 
59 /***********************************************************************
60  *
61  * Code generated by the TAYLOR translator.
62  */
63 
64 #define _N_DIM_ 9
65 
66 
67 
68 
69 
70 
71 
72 /*
73 next line defines the largest power of 2 such that 2^(LEXP2) and
74 2^(-LEXP2) do not overflow/underflow the double arithmetic of your
75 computer.
76 */
77 #define LEXP2 1023
78 
79 #define DEBUG_LEVEL 0 /* to print some internal information */
80 
81 int taylor_step_kepler_info(MY_FLOAT *ti,
82  MY_FLOAT *x,
83  int dir,
84  int step_ctl,
85  double log10abserr,
86  double log10relerr,
87  MY_FLOAT *endtime,
88  MY_FLOAT *ht,
89  int *order)
90 /*
91  * single integration step with taylor method. the parameters are:
92  *
93  * ti: on input: time of the initial condition
94  * on output: new time
95  *
96  * x: on input: initial condition
97  * on output: new condition, corresponding to the (output) time ti
98  *
99  * dir: flag to integrate forward or backwards.
100  * 1: forward
101  * -1: backwards
102  * WARNING: this flag is ignored if step_ctl (see below) is set to 0.
103  *
104  * step_ctl: flag for the step size control. the possible values are:
105  * 0: no step size control, so the step and order are provided by
106  * the user. the parameter ht is used as step, and the parameter
107  * order (see below) is used as the order.
108  * 1: standard stepsize control. it uses an approximation to the
109  * optimal order and to the radius of convergence of the series
110  * to approximate the 'optimal' step size. It tries to keep
111  * either the absolute or the relative errors below the given
112  * values. See the paper for more details.
113  * 2: as 1, but adding an extra condition on the stepsize h: the
114  * terms of the series --after being multiplied by the suitable
115  * power of h-- cannot grow.
116  * 3: user defined stepsize control. the code has to be included
117  * in the routine compute_timestep_user_defined (you can find
118  * this routine below). The user should also provide code for
119  * the selection of degree (see the function
120  * compute_order_user_defined below).
121  *
122  * log10abserr: decimal log of the absolute accuracy. the routine
123  * tries to keep either the absolute or the relative errors below
124  * the given thresholds.
125  *
126  * log10relerr: decimal log of the relative accuracy. the routine
127  * tries to keep either the absolute or the relative errors below
128  * the given thresholds.
129  *
130  * endtime: if NULL, it is ignored. if step_ctl (see above) is set
131  * to 0, it is also ignored. otherwise, if the step size is too
132  * large, it is reduced so that the new time ti is exactly endtime
133  * (in that case, the function returns 1).
134  *
135  * ht: on input: ignored/used as a time step (see parameter step_ctl)
136  * on output: time step used
137  *
138  * order: degree of the taylor expansion.
139  * input: this parameter is only used if step_ctl is 0,
140  * or if you add the code for the case step_ctl=3.
141  * its possible values are:
142  * < 2: the program will select degree 2 (if step_ctl is 0).
143  * >=2: the program will use this degree (if step_ctl is 0).
144  * ouput: degree used.
145  *
146  * return value:
147  * 0: ok.
148  * 1: ok, and ti=endtime. */
149 {
150  MY_FLOAT **taylor_coefficients_kepler_info(MY_FLOAT, MY_FLOAT*, int);
151  MY_FLOAT **taylor_coefficients_kepler_infoA(MY_FLOAT, MY_FLOAT*, int, int);
152  int compute_order_1_kepler_info(double, double, double, int*);
153  int comp_order_other_kepler_info(double, double, double);
154  double compute_stepsize_1_kepler_info(MY_FLOAT**, int, double, int);
155  double compute_stepsize_2_kepler_info(MY_FLOAT**, int, double, int);
156  double comp_stepsize_other_kepler_info(MY_FLOAT**, int, int, double, double, double);
157 
158  static MY_FLOAT **s,h,mtmp;
159  double xi,xnorm,dh;
160  int i,j,k,nt,flag_endtime,flag_err;
161  static int init=0;
162 
163  if (init == 0) /* initialization of MY_FLOAT variables */
164  {
165  init=1;
166  InitMyFloat(h);
167  InitMyFloat(mtmp);
168  }
169 /*
170  sup norm of the initial condition
171 */
172  xnorm=0;
173  if (step_ctl != 0)
174  {
175  for (i=0; i<_N_DIM_; i++)
176  {
177  MyFloatToDouble(xi,x[i]);
178  xi=fabs(xi);
179  if (xi > xnorm) xnorm=xi;
180  }
181  }
182 /*
183  we determine the degree of the taylor expansion.
184  this value will be stored in the variable nt.
185  the flag flag_err returns a value indicating if we are using an
186  absolute error estimate (value 1) or a relative one (value 2).
187 */
188  switch(step_ctl)
189  {
190  case 0: /* no step size control, fixed degree; both given by the user */
191  nt=(*order<2)? 2: *order; /* 2 is the minimum order allowed */
192  break;
193  case 1:
194  nt=compute_order_1_kepler_info(xnorm,log10abserr,log10relerr,&flag_err);
195  break;
196  case 2:
197  nt=compute_order_1_kepler_info(xnorm,log10abserr,log10relerr,&flag_err);
198  break;
199  case 3:
200  nt=comp_order_other_kepler_info(xnorm,log10abserr,log10relerr);
201  break;
202  default:
203  fprintf(stderr, "taylor_step: undefined step size control.\n");
204  fprintf(stderr, "you must choose an existing step size control\n");
205  fprintf(stderr, "method or supply a step size control procedure.\n");
206  exit(0);
207  }
208  *order=nt;
209 /*
210  computation of the jet of derivatives up to order nt
211 */
212  if(step_ctl != 0) {
213  s=taylor_coefficients_kepler_infoA(*ti,x,nt,1);
214  } else {
215  s=taylor_coefficients_kepler_info(*ti,x,nt);
216  }
217 
218 /*
219  selection of the routine to compute the time step. the value
220  step_ctl=3 has been reserved for the user, in case she/he wants to
221  code a different method.
222 */
223  switch(step_ctl)
224  {
225  case 0: /* no step size control (fixed step size, given by the user) */
226  AssignMyFloat(h,*ht);
227  break;
228  case 1:
229  dh=compute_stepsize_1_kepler_info(s,nt,xnorm,flag_err);
230  MakeMyFloatA(h, dh);
231  break;
232  case 2:
233  dh=compute_stepsize_2_kepler_info(s,nt,xnorm,flag_err);
234  MakeMyFloatA(h, dh);
235  break;
236  case 3:
237  dh=comp_stepsize_other_kepler_info(s,_N_DIM_,nt,xnorm,log10abserr,log10relerr);
238  MakeMyFloatA(h, dh);
239  break;
240  default:
241  fprintf(stderr, "taylor_step: undefined step size control.\n");
242  fprintf(stderr, "You must choose an existing step size control\n");
243  fprintf(stderr, "method or supply a step size control procedure.\n");
244  exit(0);
245  }
246 /*
247  if step_ctl != 0, we adjust the sign of the computed stepsize.
248 */
249  flag_endtime=0;
250  if (step_ctl != 0)
251  {
252  if (dir == -1) { NegateMyFloatA(mtmp,h); AssignMyFloat(h, mtmp);}
253 /*
254  we compare *ti+h with endtime. we modify h if necessary.
255 */
256  if (endtime != NULL)
257  {
258  AddMyFloatA(mtmp,h,*ti);
259  if (dir == 1) /* time goes forward */
260  {
261  if (MyFloatA_GE_B(mtmp,*endtime))
262  {
263  SubstractMyFloatA(h,*endtime,*ti);
264  flag_endtime=1;
265  }
266  }
267  else /* time goes backwards */
268  {
269  if (MyFloatA_GE_B(*endtime,mtmp))
270  {
271  SubstractMyFloatA(h,*endtime,*ti);
272  flag_endtime=1;
273  }
274  }
275  }
276  }
277 /*
278  next lines are the summation of the taylor series (horner's method)
279 */
280  j=nt-1;
281  for(i=0; i<_N_DIM_; i++) {AssignMyFloat(x[i],s[i][nt]);}
282  for(k=j; k>=0; k--)
283  {
284  for(i=0; i<_N_DIM_; i++)
285  {
286  MultiplyMyFloatA(mtmp, h, x[i]);
287  AddMyFloatA(x[i], mtmp, s[i][k]);
288  }
289  }
290 /*
291  finally, we set the values of the parameters *ht and *ti.
292 */
293  AssignMyFloat(*ht,h);
294  if (flag_endtime == 0)
295  {
296  AssignMyFloat(mtmp, *ti);
297  AddMyFloatA(*ti, mtmp, h);
298  }
299  else
300  {
301  AssignMyFloat(*ti,*endtime);
302  }
303  return(flag_endtime);
304 }
305 int compute_order_1_kepler_info(double xnorm, double log10abserr, double log10relerr, int* flag_err)
306 /*
307  * this is to determine the 'optimal' degree of the taylor expansion.
308  *
309  * parameters:
310  * xnorm: norm of the initial condition
311  * log10abserr: base-10 log of the absolute error required
312  * log10relerr: base-10 log of the relative error required
313  * flag_err: flag. returns 1 if absolute error is used
314  * returns 2 if relative error is used
315  *
316  * returned value: 'optimal' degree.
317 */
318 {
319  double log10eps,tmp;
320  int nt;
321 
322  log10eps=log10abserr;
323  *flag_err=1;
324  if (xnorm != (double)0.0)
325  {
326  tmp=log10relerr+log10(xnorm);
327  if (tmp > log10abserr) {log10eps=log10relerr; *flag_err=2;}
328  }
329 /*
330  we use 1.16 for the value 0.5*log(10)=1.151292546497...
331 */
332  nt=(int)(1.5-1.16*log10eps);
333  if (nt < 2) nt=2; /* this is the minimum order accepted */
334 
335 #if DEBUG_LEVEL > 0
336  fprintf(stderr, "taylor_step: order is %d\n",nt);
337 #endif
338 
339  return(nt);
340 }
341 double compute_stepsize_1_kepler_info(MY_FLOAT **s, int nt, double xnorm, int flag_err)
342 /*
343  * it looks for a step size for an expansion up to order nt. this
344  * function requires that nt is the value computed by
345  * compute_order_1_
346  */
347 {
348  double double_log_MyFloat_kepler_info(MY_FLOAT x);
349  static MY_FLOAT z,v1,v2;
350  static MY_FLOAT of,uf;
351  double lnv1,lnv2,r,lnro1,lnro2,lnro;
352  int i;
353  static int init=0;
354 
355  if (init == 0)
356  {
357  init=1;
358  InitMyFloat(z);
359  InitMyFloat(v1);
360  InitMyFloat(v2);
361  InitMyFloat(of);
362  InitMyFloat(uf);
363 
364  r=pow((double)2,(double)LEXP2);
365  MakeMyFloatA(of,r);
366  r=pow((double)2,(double)(-LEXP2));
367  MakeMyFloatA(uf,r);
368  }
369 /*
370  we compute the sup norm of the last two coefficients of the taylor
371  series, and we store them into v1 and v2.
372 */
373  fabsMyFloatA(v1,s[0][nt-1]);
374  fabsMyFloatA(v2,s[0][nt]);
375  for(i=1; i<_N_DIM_; i++)
376  {
377  fabsMyFloatA(z,s[i][nt-1]);
378  if (MyFloatA_GT_B(z,v1)) AssignMyFloat(v1,z);
379  fabsMyFloatA(z,s[i][nt]);
380  if (MyFloatA_GT_B(z,v2)) AssignMyFloat(v2,z);
381  }
382 /*
383  computation of the step size. we need the logs of v1 and v2, in
384  double precision (there is no need for extended precision). the idea
385  is to assign these variables to double variables and then to use the
386  standard log function. before doing that, we have to be sure that v1
387  can be assigned to a double without under or overflow. for this
388  reason we will check for this condition and, if it fails, we will
389  call an specific function for this case.
390 */
391  if (MyFloatA_LE_B(v1,of) && MyFloatA_GE_B(v1,uf))
392  {
393  MyFloatToDouble(r,v1);
394  lnv1=log(r);
395  }
396  else
397  {
398  lnv1=double_log_MyFloat_kepler_info(v1);
399  }
400  if (MyFloatA_LE_B(v2,of) && MyFloatA_GE_B(v2,uf))
401  {
402  MyFloatToDouble(r,v2);
403  lnv2=log(r);
404  }
405  else
406  {
407  lnv2=double_log_MyFloat_kepler_info(v2);
408  }
409 /*
410  if flag_err == 2, this means that we are using a relative error control.
411 */
412  if (flag_err == 2)
413  {
414  r = -log10(xnorm);
415  lnv1 += r;
416  lnv2 += r;
417  }
418  lnro1= -lnv1/(nt-1);
419  lnro2= -lnv2/nt;
420  lnro=(lnro1 < lnro2)? lnro1: lnro2;
421 
422  r=exp(lnro-2-0.7/(nt-1)); /* exp(-0.7/(nt-1)) is a security factor */
423  return(r);
424 }
425 double compute_stepsize_2_kepler_info(MY_FLOAT **s, int nt, double xnorm, int flag_err)
426 /*
427  * it looks for a step size for an expansion up to order nt. this
428  * function requires that nt is the value computed by
429  * compute_order_1_. it also tries to reduce cancellations of big
430  * terms in the summation of the taylor series.
431  */
432 {
433  double compute_stepsize_1_kepler_info(MY_FLOAT**, int, double, int);
434  static MY_FLOAT h,hj,r,z,a,normj;
435  double c,rtmp,dh;
436  int i,j;
437  static int init=0;
438 
439  if (init == 0)
440  {
441  init=1;
442  InitMyFloat(h);
443  InitMyFloat(hj);
444  InitMyFloat(r);
445  InitMyFloat(z);
446  InitMyFloat(a);
447  InitMyFloat(normj);
448  }
449 /*
450  we compute the step size according to the first algorithm
451 */
452  dh=compute_stepsize_1_kepler_info(s,nt,xnorm,flag_err);
453  MakeMyFloatA(h,dh);
454 /*
455  next lines select a value (z), that will be used to control the size
456  of the terms of the Taylor series.
457 */
458  if (flag_err == 1) {
459  MakeMyFloatA(z, 1.0);
460  } else if (flag_err == 2) {
461  MakeMyFloatA(z,xnorm);
462  } else
463  {
464  printf("compute_stepsize_2 internal error. flag_err: %d\n",flag_err);
465  exit(1);
466  }
467 /*
468  next loop checks if the sup norm of the terms in the Taylor series are
469  lower than z. if a term is greater than z, the step size h is reduced.
470 */
471  MakeMyFloatA(hj,(double)1.0);
472 
473  for(j=1; j<=nt; j++)
474  {
475  MultiplyMyFloatA(r,h,hj);
476  AssignMyFloat(hj,r);
477 
478  MakeMyFloatC(normj,"0", (double)0);
479  for (i=0; i<_N_DIM_; i++)
480  {
481  fabsMyFloatA(a,s[i][j]);
482  if (MyFloatA_GT_B(a,normj)) AssignMyFloat(normj,a);
483  }
484 
485  MultiplyMyFloatA(r,normj,hj);
486  if (MyFloatA_LE_B(r,z)) continue;
487 /*
488  we reduce h (and hj)
489 */
490  DivideMyFloatA(hj,z,normj);
491 
492  DivideMyFloatA(a,r,z);
493  MyFloatToDouble(c,a);
494  c=pow(c,(double)1.e0/(double)j);
495  MakeMyFloatA(a,c);
496  DivideMyFloatA(r,h,a);
497  AssignMyFloat(h,r);
498 
499 #if DEBUG_LEVEL > 1
500  fprintf(stderr, "order %2d. reducing h from %14.6e to %14.6e\n",j,c*h,h);
501 #endif
502  }
503 
504  MyFloatToDouble(rtmp,h);
505  return(rtmp);
506 }
507 
508 double double_log_MyFloat_kepler_info(MY_FLOAT x)
509 /*
510  * natural log, in double precision, of a MY_FLOAT positive number.
511  */
512 {
513  static MY_FLOAT a,tmp;
514  static MY_FLOAT z,of,uf;
515  double b,lx;
516  int k;
517  static int init=0;
518 
519  if (init == 0)
520  {
521  init=1;
522  InitMyFloat(a);
523  InitMyFloat(z);
524  InitMyFloat(of);
525  InitMyFloat(uf);
526  InitMyFloat(tmp);
527 
528  b=0;
529  MakeMyFloatA(z,b);
530  b=pow((double)2,(double)LEXP2);
531  MakeMyFloatA(of,b);
532  b=pow((double)2,(double)(-LEXP2));
533  MakeMyFloatA(uf,b);
534  }
535 
536  if (MyFloatA_EQ_B(x,z))
537  {
538  puts("double_log_MyFloat error: zero argument");
539  puts("(this is because one of the last two terms of your taylor");
540  puts(" expansion is exactly zero)");
541  exit(1);
542  }
543 
544  AssignMyFloat(a,x);
545 
546  k=0;
547  while(MyFloatA_LT_B(a,uf))
548  {
549  ++k;
550  if(k>3000){fprintf(stderr,"double_log_MyFloat overflow: %d\n", k); exit(1);}
551  MultiplyMyFloatA(tmp,a,of);
552  AssignMyFloat(a,tmp);
553  }
554  while(MyFloatA_GT_B(a,of))
555  {
556  --k;
557  if(k<-3000){fprintf(stderr,"double_log_MyFloat underflow: %d\n", k); exit(1);}
558  MultiplyMyFloatA(tmp,a,uf);
559  AssignMyFloat(a,tmp);
560  }
561 
562  MyFloatToDouble(b,a);
563 /*
564  lx stands for log(x)
565 */
566  lx=log(b)-(LEXP2*0.69314718055994530942)*k;
567 
568  return(lx);
569 }
570 
571 
572 int comp_order_other_kepler_info(double lnxnorm, double log10abserr, double log10relerr){
573  puts("---");
574  puts("compute_order_user_defined:");
575  puts("you have to code this routine");
576  puts("or select a different value for the step_ctl parameter");
577  puts("---");
578  exit(1);
579 
580  return(0);
581 }
582 double comp_stepsize_other_kepler_info(MY_FLOAT **s, int nd, int nt, double xnorm, double log10abserr, double log10relerr) {
583 
584  puts("---");
585  puts("compute_timestep_user_defined:");
586  puts("you have to code this routine");
587  puts("or select a different value for the step_ctl parameter");
588  puts("---");
589  exit(1);
590  return((double)0.00001);
591 }
592 /***********************************************************************
593  *
594  * Procedure generated by the TAYLOR translator. Do not edit!
595  *
596  * It needs the header file 'taylor.h' to compile.
597  * Run taylor with the -header -o taylor.h option to generate a sample 'taylor.h'
598 
599  * Translation info is at the end of this file.
600  * Version 1.4.4, Apr 22, 2008
601  ***********************************************************************/
602 
603 #include <stdio.h>
604 #include <stdlib.h>
605 MY_FLOAT **taylor_coefficients_kepler_infoA(MY_FLOAT t, MY_FLOAT *x, int order, int rflag)
606 {
607  /* input:
608  t: current value of the time variable
609  x: array represent values of the state variables
610  order: order of the taylor coefficients sought
611  rflag: recompute flag. If you call this routine with one order
612  first, but then decided that you need a higher order of the
613  taylor polynomial. You can pass 0 to rflag. This routine
614  will try to use the values already computed. Provided that
615  both x and t have not been changed, and you did not modify
616  the jet derivatives from the previous call.
617  Return Value:
618  Two D Array, rows are the taylor coefficients of the
619  state variables
620 
621  */
622 
623  static int _jz_ivars[2];
624  static MY_FLOAT *_jz_jet[28], *_jz_save = NULL, *_jz_oneOverN=NULL,*_jz_theNs=NULL;
625  static MY_FLOAT _jz_tvar1, _jz_tvar2, _jz_tvar3, _jz_tvar4; /* tmp vars */
626  static MY_FLOAT _jz_uvar1, _jz_uvar2; /* tmp vars */
627  static MY_FLOAT _jz_svar1, _jz_svar2, _jz_svar3, _jz_svar4, _jz_svar5; /* tmp vars */
628  static MY_FLOAT _jz_wvar3, _jz_wvar4; /* tmp vars */
629  static MY_FLOAT _jz_zvar1, _jz_zvar2; /* tmp vars */
630  static MY_FLOAT _jz_MyFloatZERO;
631  static int _jz_maxOrderUsed = -1;
632  static int _jz_lastOrder = 0, _jz_initialized=0, _jz_ginitialized=0;
633  int _jz_i, _jz_j, _jz_k, _jz_l, _jz_m, _jz_n, _jz_oorder ;
634  /* allocating memory if needed */
635  if(_jz_maxOrderUsed < order ) {
636  if(_jz_ginitialized == 0) {
637  InitMyFloat(_jz_tvar1); InitMyFloat(_jz_tvar2);InitMyFloat(_jz_tvar3);InitMyFloat(_jz_tvar4);
638  InitMyFloat(_jz_svar1); InitMyFloat(_jz_svar2);InitMyFloat(_jz_svar3);InitMyFloat(_jz_svar4);
639  InitMyFloat(_jz_svar5); InitMyFloat(_jz_zvar1);InitMyFloat(_jz_zvar2);
640  InitMyFloat(_jz_uvar1); InitMyFloat(_jz_uvar2);
641  InitMyFloat(_jz_wvar3);InitMyFloat(_jz_wvar4);
642  InitMyFloat(_jz_MyFloatZERO);
643  MakeMyFloatC(_jz_MyFloatZERO, "0", (double)0);
644  }
645  if(rflag > 0) rflag = 0; /* have to recompute everything */
646  _jz_oorder=_jz_maxOrderUsed;
647  _jz_maxOrderUsed = order;
648  if(_jz_ginitialized) {
649  for(_jz_i=0; _jz_i< _jz_oorder+1; _jz_i++) {ClearMyFloat(_jz_oneOverN[_jz_i]); ClearMyFloat(_jz_theNs[_jz_i]);} free(_jz_oneOverN); free(_jz_theNs);
650  }
651  _jz_theNs = (MY_FLOAT *)malloc((order+1) * sizeof(MY_FLOAT));
652  _jz_oneOverN = (MY_FLOAT *)malloc((order+1) * sizeof(MY_FLOAT));
653  for(_jz_i=0; _jz_i<order+1; _jz_i++) {InitMyFloat(_jz_oneOverN[_jz_i]);InitMyFloat(_jz_theNs[_jz_i]);}
654  MakeMyFloatC(_jz_theNs[0],"0.0", (double)0.0);
655  MakeMyFloatC(_jz_uvar1,"1.0", (double)1.0);
656  for(_jz_i = 1; _jz_i <= order; _jz_i++) {
657  AssignMyFloat(_jz_tvar2, _jz_theNs[_jz_i-1]);
658  AddMyFloatA(_jz_theNs[_jz_i], _jz_tvar2, _jz_uvar1);
659  }
660  AssignMyFloat(_jz_oneOverN[0],_jz_uvar1);
661  AssignMyFloat(_jz_oneOverN[1],_jz_uvar1);
662  for(_jz_i = 2; _jz_i <= order; _jz_i++) {
663  DivideMyFloatA(_jz_oneOverN[_jz_i], _jz_uvar1,_jz_theNs[_jz_i]);
664  }
665  if(_jz_ginitialized) {
666  for(_jz_i=0; _jz_i<(_jz_oorder+1)*(28); _jz_i++) { ClearMyFloat(_jz_save[_jz_i]);} free(_jz_save);
667  }
668  _jz_save = (MY_FLOAT *)malloc((order+1)* 28 *sizeof(MY_FLOAT));
669  for(_jz_i=0; _jz_i<(order+1)*(28); _jz_i++) { InitMyFloat(_jz_save[_jz_i]);}
670  for(_jz_j = 0, _jz_k = 0; _jz_j < 28 ; _jz_j++, _jz_k += order+1) { _jz_jet[_jz_j] =& (_jz_save[_jz_k]); }
671 
672  /* True constants, initialized only once. */
673  /* const: i_022=2 */
674  _jz_ivars[0]=2;
675  /* const: i_026=3 */
676  _jz_ivars[1]=3;
677  /* const: v_029=0 */
678  MakeMyFloatC(_jz_jet[25][0],"0",(double)0);
679  /* const: v_030=0 */
680  MakeMyFloatC(_jz_jet[26][0],"0",(double)0);
681  /* const: v_031=0 */
682  MakeMyFloatC(_jz_jet[27][0],"0",(double)0);
683  }
684 
685  if(rflag) {
686  if(rflag < 0 ) return(NULL);
687  for(_jz_i = 0; rflag != 0 && _jz_i < 9; _jz_i++) {
688  if(MyFloatA_NEQ_B(_jz_jet[_jz_i][0], x[_jz_i])) rflag = 0;
689  }
690  }
691 
692  if(rflag == 0) {
693  /* initialize all constant vars and state variables */
694  _jz_lastOrder = 1;
695  AssignMyFloat(_jz_jet[0][0], x[0]);
696  AssignMyFloat(_jz_jet[1][0], x[1]);
697  AssignMyFloat(_jz_jet[2][0], x[2]);
698  AssignMyFloat(_jz_jet[3][0], x[3]);
699  AssignMyFloat(_jz_jet[4][0], x[4]);
700  AssignMyFloat(_jz_jet[5][0], x[5]);
701  AssignMyFloat(_jz_jet[6][0], x[6]);
702  AssignMyFloat(_jz_jet[7][0], x[7]);
703  AssignMyFloat(_jz_jet[8][0], x[8]);
704  /* negate: v_038=(-v_011) */
705  NegateMyFloatA(_jz_jet[9][0],_jz_jet[0][0]);
706  /* exponentiation: v_032=(v_011^i_022) */
707  /* integer exponent or half integer */
708  AssignMyFloat(_jz_svar5,_jz_jet[0][0]);
709  { int n=2, m, mn=0;
710  switch(n) {
711  case 0: AssignMyFloat(_jz_jet[10][0], _jz_oneOverN[0]); break;
712  case 1: AssignMyFloat(_jz_jet[10][0], _jz_svar5); break;
713  case 2: AssignMyFloat(_jz_svar1, _jz_svar5); MultiplyMyFloatA(_jz_jet[10][0],_jz_svar1,_jz_svar5); break;
714  case 3: AssignMyFloat(_jz_svar1, _jz_svar5); MultiplyMyFloatA(_jz_svar2,_jz_svar1,_jz_svar5);
715  MultiplyMyFloatA(_jz_jet[10][0],_jz_svar1,_jz_svar2); break;
716  default:
717  AssignMyFloat(_jz_svar1, _jz_oneOverN[0]); AssignMyFloat(_jz_svar2, _jz_svar5);
718  while(mn==0) {
719  m=n; n /=2; if(n+n != m) {
720  AssignMyFloat(_jz_svar3, _jz_svar1); MultiplyMyFloatA(_jz_svar1, _jz_svar3, _jz_svar2);
721  if(n==0){ mn=1; AssignMyFloat(_jz_jet[10][0],_jz_svar1);}
722  }
723  if(mn==0) {AssignMyFloat(_jz_svar3, _jz_svar2);AssignMyFloat(_jz_svar4, _jz_svar2);
724  MultiplyMyFloatA(_jz_svar2, _jz_svar3,_jz_svar4);}
725  }
726  break;
727  }
728  }
729  /* exponentiation: v_033=(v_012^i_022) */
730  /* integer exponent or half integer */
731  AssignMyFloat(_jz_svar5,_jz_jet[1][0]);
732  { int n=2, m, mn=0;
733  switch(n) {
734  case 0: AssignMyFloat(_jz_jet[11][0], _jz_oneOverN[0]); break;
735  case 1: AssignMyFloat(_jz_jet[11][0], _jz_svar5); break;
736  case 2: AssignMyFloat(_jz_svar1, _jz_svar5); MultiplyMyFloatA(_jz_jet[11][0],_jz_svar1,_jz_svar5); break;
737  case 3: AssignMyFloat(_jz_svar1, _jz_svar5); MultiplyMyFloatA(_jz_svar2,_jz_svar1,_jz_svar5);
738  MultiplyMyFloatA(_jz_jet[11][0],_jz_svar1,_jz_svar2); break;
739  default:
740  AssignMyFloat(_jz_svar1, _jz_oneOverN[0]); AssignMyFloat(_jz_svar2, _jz_svar5);
741  while(mn==0) {
742  m=n; n /=2; if(n+n != m) {
743  AssignMyFloat(_jz_svar3, _jz_svar1); MultiplyMyFloatA(_jz_svar1, _jz_svar3, _jz_svar2);
744  if(n==0){ mn=1; AssignMyFloat(_jz_jet[11][0],_jz_svar1);}
745  }
746  if(mn==0) {AssignMyFloat(_jz_svar3, _jz_svar2);AssignMyFloat(_jz_svar4, _jz_svar2);
747  MultiplyMyFloatA(_jz_svar2, _jz_svar3,_jz_svar4);}
748  }
749  break;
750  }
751  }
752  /* plus: v_034=(v_032+v_033) */
753  AddMyFloatA(_jz_jet[12][0], _jz_jet[10][0], _jz_jet[11][0]);
754  /* exponentiation: v_035=(v_013^i_022) */
755  /* integer exponent or half integer */
756  AssignMyFloat(_jz_svar5,_jz_jet[2][0]);
757  { int n=2, m, mn=0;
758  switch(n) {
759  case 0: AssignMyFloat(_jz_jet[13][0], _jz_oneOverN[0]); break;
760  case 1: AssignMyFloat(_jz_jet[13][0], _jz_svar5); break;
761  case 2: AssignMyFloat(_jz_svar1, _jz_svar5); MultiplyMyFloatA(_jz_jet[13][0],_jz_svar1,_jz_svar5); break;
762  case 3: AssignMyFloat(_jz_svar1, _jz_svar5); MultiplyMyFloatA(_jz_svar2,_jz_svar1,_jz_svar5);
763  MultiplyMyFloatA(_jz_jet[13][0],_jz_svar1,_jz_svar2); break;
764  default:
765  AssignMyFloat(_jz_svar1, _jz_oneOverN[0]); AssignMyFloat(_jz_svar2, _jz_svar5);
766  while(mn==0) {
767  m=n; n /=2; if(n+n != m) {
768  AssignMyFloat(_jz_svar3, _jz_svar1); MultiplyMyFloatA(_jz_svar1, _jz_svar3, _jz_svar2);
769  if(n==0){ mn=1; AssignMyFloat(_jz_jet[13][0],_jz_svar1);}
770  }
771  if(mn==0) {AssignMyFloat(_jz_svar3, _jz_svar2);AssignMyFloat(_jz_svar4, _jz_svar2);
772  MultiplyMyFloatA(_jz_svar2, _jz_svar3,_jz_svar4);}
773  }
774  break;
775  }
776  }
777  /* plus: v_036=(v_034+v_035) */
778  AddMyFloatA(_jz_jet[14][0], _jz_jet[12][0], _jz_jet[13][0]);
779  /* call: v_037=sqrt(v_036) */
780  sqrtMyFloatA(_jz_jet[15][0], _jz_jet[14][0]);
781  /* exponentiation: v_039=(v_037^i_026) */
782  /* integer exponent or half integer */
783  AssignMyFloat(_jz_svar5,_jz_jet[15][0]);
784  { int n=3, m, mn=0;
785  switch(n) {
786  case 0: AssignMyFloat(_jz_jet[16][0], _jz_oneOverN[0]); break;
787  case 1: AssignMyFloat(_jz_jet[16][0], _jz_svar5); break;
788  case 2: AssignMyFloat(_jz_svar1, _jz_svar5); MultiplyMyFloatA(_jz_jet[16][0],_jz_svar1,_jz_svar5); break;
789  case 3: AssignMyFloat(_jz_svar1, _jz_svar5); MultiplyMyFloatA(_jz_svar2,_jz_svar1,_jz_svar5);
790  MultiplyMyFloatA(_jz_jet[16][0],_jz_svar1,_jz_svar2); break;
791  default:
792  AssignMyFloat(_jz_svar1, _jz_oneOverN[0]); AssignMyFloat(_jz_svar2, _jz_svar5);
793  while(mn==0) {
794  m=n; n /=2; if(n+n != m) {
795  AssignMyFloat(_jz_svar3, _jz_svar1); MultiplyMyFloatA(_jz_svar1, _jz_svar3, _jz_svar2);
796  if(n==0){ mn=1; AssignMyFloat(_jz_jet[16][0],_jz_svar1);}
797  }
798  if(mn==0) {AssignMyFloat(_jz_svar3, _jz_svar2);AssignMyFloat(_jz_svar4, _jz_svar2);
799  MultiplyMyFloatA(_jz_svar2, _jz_svar3,_jz_svar4);}
800  }
801  break;
802  }
803  }
804  /* div: v_040=(v_038/v_039) */
805  DivideMyFloatA(_jz_jet[17][0], _jz_jet[9][0], _jz_jet[16][0]);
806  /* plus: v_041=(v_017+v_040) */
807  AddMyFloatA(_jz_jet[18][0], _jz_jet[6][0], _jz_jet[17][0]);
808  /* negate: v_042=(-v_012) */
809  NegateMyFloatA(_jz_jet[19][0],_jz_jet[1][0]);
810  /* div: v_044=(v_042/v_039) */
811  DivideMyFloatA(_jz_jet[20][0], _jz_jet[19][0], _jz_jet[16][0]);
812  /* plus: v_045=(v_018+v_044) */
813  AddMyFloatA(_jz_jet[21][0], _jz_jet[7][0], _jz_jet[20][0]);
814  /* negate: v_046=(-v_013) */
815  NegateMyFloatA(_jz_jet[22][0],_jz_jet[2][0]);
816  /* div: v_048=(v_046/v_039) */
817  DivideMyFloatA(_jz_jet[23][0], _jz_jet[22][0], _jz_jet[16][0]);
818  /* plus: v_049=(v_019+v_048) */
819  AddMyFloatA(_jz_jet[24][0], _jz_jet[8][0], _jz_jet[23][0]);
820 
821  /* the first derivative of state variables */
822  /* state variable 0: */
823  AssignMyFloat(_jz_jet[0][1], _jz_jet[3][0]);
824  /* state variable 1: */
825  AssignMyFloat(_jz_jet[1][1], _jz_jet[4][0]);
826  /* state variable 2: */
827  AssignMyFloat(_jz_jet[2][1], _jz_jet[5][0]);
828  /* state variable 3: */
829  AssignMyFloat(_jz_jet[3][1], _jz_jet[18][0]);
830  /* state variable 4: */
831  AssignMyFloat(_jz_jet[4][1], _jz_jet[21][0]);
832  /* state variable 5: */
833  AssignMyFloat(_jz_jet[5][1], _jz_jet[24][0]);
834  /* state variable 6: */
835  AssignMyFloat(_jz_jet[6][1], _jz_jet[25][0]);
836  /* state variable 7: */
837  AssignMyFloat(_jz_jet[7][1], _jz_jet[26][0]);
838  /* state variable 8: */
839  AssignMyFloat(_jz_jet[8][1], _jz_jet[27][0]);
840  }
841 
842  /* compute the kth order derivatives of all vars */
843  for(_jz_k = _jz_lastOrder; _jz_k < order; _jz_k++) {
844  /* derivative for tmp variables */
845  /* negation: v_038=(-v_011) */
846  NegateMyFloatA(_jz_jet[9][_jz_k], _jz_jet[0][_jz_k]);
847  /* exponentiation: v_032=(v_011^i_022) */
848  { /* exponentiation */
849  /* expr^2 */
850  static MY_FLOAT tmp1, tmp2, tmp;
851  int parity=(_jz_k&1), half=(_jz_k+1)>>1;
852  if(_jz_initialized==0) { InitMyFloat(tmp1);InitMyFloat(tmp2); InitMyFloat(tmp);}
853  AssignMyFloat(tmp, _jz_MyFloatZERO);
854  for(_jz_l=0; _jz_l<half; _jz_l++) {
855  MultiplyMyFloatA(tmp1, _jz_jet[0][_jz_l], _jz_jet[0][_jz_k-_jz_l]);
856  AssignMyFloat(tmp2, tmp);
857  AddMyFloatA(tmp, tmp2, tmp1);
858  }
859  AssignMyFloat(tmp2, tmp);
860  AddMyFloatA(tmp1, tmp2, tmp);
861  if(parity==0) {
862  MultiplyMyFloatA(tmp2, _jz_jet[0][half], _jz_jet[0][half]);
863  AddMyFloatA(_jz_jet[10][_jz_k], tmp2, tmp1);
864  } else {
865  AssignMyFloat(_jz_jet[10][_jz_k], tmp1);
866  }
867  }
868  /* exponentiation: v_033=(v_012^i_022) */
869  { /* exponentiation */
870  /* expr^2 */
871  static MY_FLOAT tmp1, tmp2, tmp;
872  int parity=(_jz_k&1), half=(_jz_k+1)>>1;
873  if(_jz_initialized==0) { InitMyFloat(tmp1);InitMyFloat(tmp2); InitMyFloat(tmp);}
874  AssignMyFloat(tmp, _jz_MyFloatZERO);
875  for(_jz_l=0; _jz_l<half; _jz_l++) {
876  MultiplyMyFloatA(tmp1, _jz_jet[1][_jz_l], _jz_jet[1][_jz_k-_jz_l]);
877  AssignMyFloat(tmp2, tmp);
878  AddMyFloatA(tmp, tmp2, tmp1);
879  }
880  AssignMyFloat(tmp2, tmp);
881  AddMyFloatA(tmp1, tmp2, tmp);
882  if(parity==0) {
883  MultiplyMyFloatA(tmp2, _jz_jet[1][half], _jz_jet[1][half]);
884  AddMyFloatA(_jz_jet[11][_jz_k], tmp2, tmp1);
885  } else {
886  AssignMyFloat(_jz_jet[11][_jz_k], tmp1);
887  }
888  }
889  /* plus: v_034=(v_032+v_033) */
890  AddMyFloatA(_jz_jet[12][_jz_k], _jz_jet[10][_jz_k],_jz_jet[11][_jz_k]);
891  /* exponentiation: v_035=(v_013^i_022) */
892  { /* exponentiation */
893  /* expr^2 */
894  static MY_FLOAT tmp1, tmp2, tmp;
895  int parity=(_jz_k&1), half=(_jz_k+1)>>1;
896  if(_jz_initialized==0) { InitMyFloat(tmp1);InitMyFloat(tmp2); InitMyFloat(tmp);}
897  AssignMyFloat(tmp, _jz_MyFloatZERO);
898  for(_jz_l=0; _jz_l<half; _jz_l++) {
899  MultiplyMyFloatA(tmp1, _jz_jet[2][_jz_l], _jz_jet[2][_jz_k-_jz_l]);
900  AssignMyFloat(tmp2, tmp);
901  AddMyFloatA(tmp, tmp2, tmp1);
902  }
903  AssignMyFloat(tmp2, tmp);
904  AddMyFloatA(tmp1, tmp2, tmp);
905  if(parity==0) {
906  MultiplyMyFloatA(tmp2, _jz_jet[2][half], _jz_jet[2][half]);
907  AddMyFloatA(_jz_jet[13][_jz_k], tmp2, tmp1);
908  } else {
909  AssignMyFloat(_jz_jet[13][_jz_k], tmp1);
910  }
911  }
912  /* plus: v_036=(v_034+v_035) */
913  AddMyFloatA(_jz_jet[14][_jz_k], _jz_jet[12][_jz_k],_jz_jet[13][_jz_k]);
914  /* call: v_037=sqrt(v_036) */
915  { /* call sqrt */
916  static MY_FLOAT tmp1, tmp2, tmp3, tmp;
917  if(_jz_initialized==0) {
918  InitMyFloat(tmp);InitMyFloat(tmp1);InitMyFloat(tmp2);
919  InitMyFloat(tmp3);
920  }
921  AssignMyFloat(tmp, _jz_jet[14][_jz_k]);
922  for(_jz_l=1; _jz_l< _jz_k; _jz_l++) {
923  MultiplyMyFloatA(tmp1, _jz_jet[15][_jz_k-_jz_l],_jz_jet[15][_jz_l]);
924  SubstractMyFloatA(tmp3, tmp, tmp1);
925  AssignMyFloat(tmp, tmp3);
926  }
927  AssignMyFloat(tmp1, _jz_jet[15][0]);
928  AddMyFloatA(tmp3, tmp1, tmp1);
929  DivideMyFloatA(_jz_jet[15][_jz_k], tmp, tmp3);
930  }
931  /* exponentiation: v_039=(v_037^i_026) */
932  { /* exponentiation */
933  /* expr^(3)/ */
934  int ppk=(3)*_jz_k, qqk=_jz_k, pq=(3+1);
935  static MY_FLOAT tmp1, tmp2, tmp3, tmpC, tmp;
936  if(_jz_initialized==0) {
937  InitMyFloat(tmp1);InitMyFloat(tmp2); InitMyFloat(tmp3);
938  InitMyFloat(tmpC);InitMyFloat(tmp);
939  }
940  AssignMyFloat(tmp, _jz_MyFloatZERO);
941  for(_jz_l=0; _jz_l<_jz_k; _jz_l++) {
942  MakeMyFloatA(tmpC, ppk);
943  ppk -= pq ;
944  MultiplyMyFloatA(tmp1, _jz_jet[16][_jz_l],_jz_jet[15][_jz_k-_jz_l]);
945  MultiplyMyFloatA(tmp2, tmpC, tmp1);
946  AddMyFloatA(tmp1, tmp, tmp2);
947  AssignMyFloat(tmp, tmp1);
948  }
949  MakeMyFloatA(tmp3,qqk);
950  MultiplyMyFloatA(tmp1, _jz_jet[15][0], tmp3);
951  DivideMyFloatA(_jz_jet[16][_jz_k], tmp, tmp1);
952  }
953  /* div: v_040=(v_038/v_039) */
954  { /* division */
955  static MY_FLOAT tmp1, tmp2, tmp;
956  if(_jz_initialized==0) { InitMyFloat(tmp1);InitMyFloat(tmp2); InitMyFloat(tmp);}
957  AssignMyFloat(tmp, _jz_MyFloatZERO);
958  for(_jz_l=1; _jz_l<=_jz_k; _jz_l++) {
959  MultiplyMyFloatA(tmp1, _jz_jet[16][_jz_l],_jz_jet[17][_jz_k-_jz_l]);
960  AssignMyFloat(tmp2, tmp);
961  AddMyFloatA(tmp, tmp2, tmp1);
962  }
963  AssignMyFloat(tmp2, tmp);
964  SubstractMyFloatA(tmp, _jz_jet[9][_jz_k], tmp2);
965  DivideMyFloatA(_jz_jet[17][_jz_k], tmp, _jz_jet[16][0]);
966  }
967  /* plus: v_041=(v_017+v_040) */
968  AddMyFloatA(_jz_jet[18][_jz_k], _jz_jet[6][_jz_k],_jz_jet[17][_jz_k]);
969  /* negation: v_042=(-v_012) */
970  NegateMyFloatA(_jz_jet[19][_jz_k], _jz_jet[1][_jz_k]);
971  /* div: v_044=(v_042/v_039) */
972  { /* division */
973  static MY_FLOAT tmp1, tmp2, tmp;
974  if(_jz_initialized==0) { InitMyFloat(tmp1);InitMyFloat(tmp2); InitMyFloat(tmp);}
975  AssignMyFloat(tmp, _jz_MyFloatZERO);
976  for(_jz_l=1; _jz_l<=_jz_k; _jz_l++) {
977  MultiplyMyFloatA(tmp1, _jz_jet[16][_jz_l],_jz_jet[20][_jz_k-_jz_l]);
978  AssignMyFloat(tmp2, tmp);
979  AddMyFloatA(tmp, tmp2, tmp1);
980  }
981  AssignMyFloat(tmp2, tmp);
982  SubstractMyFloatA(tmp, _jz_jet[19][_jz_k], tmp2);
983  DivideMyFloatA(_jz_jet[20][_jz_k], tmp, _jz_jet[16][0]);
984  }
985  /* plus: v_045=(v_018+v_044) */
986  AddMyFloatA(_jz_jet[21][_jz_k], _jz_jet[7][_jz_k],_jz_jet[20][_jz_k]);
987  /* negation: v_046=(-v_013) */
988  NegateMyFloatA(_jz_jet[22][_jz_k], _jz_jet[2][_jz_k]);
989  /* div: v_048=(v_046/v_039) */
990  { /* division */
991  static MY_FLOAT tmp1, tmp2, tmp;
992  if(_jz_initialized==0) { InitMyFloat(tmp1);InitMyFloat(tmp2); InitMyFloat(tmp);}
993  AssignMyFloat(tmp, _jz_MyFloatZERO);
994  for(_jz_l=1; _jz_l<=_jz_k; _jz_l++) {
995  MultiplyMyFloatA(tmp1, _jz_jet[16][_jz_l],_jz_jet[23][_jz_k-_jz_l]);
996  AssignMyFloat(tmp2, tmp);
997  AddMyFloatA(tmp, tmp2, tmp1);
998  }
999  AssignMyFloat(tmp2, tmp);
1000  SubstractMyFloatA(tmp, _jz_jet[22][_jz_k], tmp2);
1001  DivideMyFloatA(_jz_jet[23][_jz_k], tmp, _jz_jet[16][0]);
1002  }
1003  /* plus: v_049=(v_019+v_048) */
1004  AddMyFloatA(_jz_jet[24][_jz_k], _jz_jet[8][_jz_k],_jz_jet[23][_jz_k]);
1005  /* constants: v_029=0 ! */
1006  AssignMyFloat(_jz_jet[25][_jz_k], _jz_MyFloatZERO);
1007  /* constants: v_030=0 ! */
1008  AssignMyFloat(_jz_jet[26][_jz_k], _jz_MyFloatZERO);
1009  /* constants: v_031=0 ! */
1010  AssignMyFloat(_jz_jet[27][_jz_k], _jz_MyFloatZERO);
1011  /* derivative of state variables */
1012  _jz_m = _jz_k+1;
1013  /* state variable 0: */
1014  DivideMyFloatByInt(_jz_jet[0][_jz_m], _jz_jet[3][_jz_k], _jz_m);
1015  /* state variable 1: */
1016  DivideMyFloatByInt(_jz_jet[1][_jz_m], _jz_jet[4][_jz_k], _jz_m);
1017  /* state variable 2: */
1018  DivideMyFloatByInt(_jz_jet[2][_jz_m], _jz_jet[5][_jz_k], _jz_m);
1019  /* state variable 3: */
1020  DivideMyFloatByInt(_jz_jet[3][_jz_m], _jz_jet[18][_jz_k], _jz_m);
1021  /* state variable 4: */
1022  DivideMyFloatByInt(_jz_jet[4][_jz_m], _jz_jet[21][_jz_k], _jz_m);
1023  /* state variable 5: */
1024  DivideMyFloatByInt(_jz_jet[5][_jz_m], _jz_jet[24][_jz_k], _jz_m);
1025  /* state variable 6: */
1026  DivideMyFloatByInt(_jz_jet[6][_jz_m], _jz_jet[25][_jz_k], _jz_m);
1027  /* state variable 7: */
1028  DivideMyFloatByInt(_jz_jet[7][_jz_m], _jz_jet[26][_jz_k], _jz_m);
1029  /* state variable 8: */
1030  DivideMyFloatByInt(_jz_jet[8][_jz_m], _jz_jet[27][_jz_k], _jz_m);
1031  _jz_initialized=1;
1032  }
1033  _jz_lastOrder = order;
1034  _jz_ginitialized=1;
1035  return(_jz_jet);
1036 }
1037 MY_FLOAT **taylor_coefficients_kepler_info(MY_FLOAT t, MY_FLOAT *x, int order)
1038 {
1039  return(taylor_coefficients_kepler_infoA(t,x,order,0));
1040 }
1041 
1042 /******************** Translation Info *****************************/
1043 /*
1044 
1045 
1046 ===================================================================================
1047 ======= ======
1048 ======= Final Variable List ======
1049  (28 + 0) vars, (0 + 0) cvars and (2 + 0) ivars
1050 ======= ======
1051 ===================================================================================
1052  v_011 (state variable)
1053  v_012 (state variable)
1054  v_013 (state variable)
1055  v_014 (state variable)
1056  v_015 (state variable)
1057  v_016 (state variable)
1058  v_017 (state variable)
1059  v_018 (state variable)
1060  v_019 (state variable)
1061  v_038 = (-v_011) (9 0)
1062  i_022 = 2 (0 0) (a number)
1063  v_032 = (v_011^i_022) (10 0)
1064  v_033 = (v_012^i_022) (11 0)
1065  v_034 = (v_032+v_033) (12 0)
1066  v_035 = (v_013^i_022) (13 0)
1067  v_036 = (v_034+v_035) (14 0)
1068  v_037 = sqrt(v_036) (15 0)
1069  i_026 = 3 (1 0) (a number)
1070  v_039 = (v_037^i_026) (16 0)
1071  v_040 = (v_038/v_039) (17 0)
1072  v_041 = (v_017+v_040) (18 0)
1073  v_042 = (-v_012) (19 0)
1074  v_044 = (v_042/v_039) (20 0)
1075  v_045 = (v_018+v_044) (21 0)
1076  v_046 = (-v_013) (22 0)
1077  v_048 = (v_046/v_039) (23 0)
1078  v_049 = (v_019+v_048) (24 0)
1079  v_029 = 0 (25 0)
1080  v_030 = 0 (26 0)
1081  v_031 = 0 (27 0)
1082 ===================================================================
1083 ========= ========
1084 ========= Differential Equations ========
1085 ========= ========
1086 ===================================================================
1087 
1088  v_011'=v_014
1089  v_012'=v_015
1090  v_013'=v_016
1091  v_014'=v_041
1092  v_015'=v_045
1093  v_016'=v_049
1094  v_017'=v_029
1095  v_018'=v_030
1096  v_019'=v_031
1097 */
1098 /*************** END END END ***************************************/
pagmo::util::neighbourhood::order
std::vector< population::size_type > order(const std::vector< T > &values)
Sort according the the values in the values vector but return the permutation.
Definition: neighbourhood.h:55
Generated by   doxygen 1.8.9.1