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