117 void write_formatted(
const double&, ostream& ) ;
118 void write_formatted(
const Tbl&, ostream& ) ;
126 template<
typename TyT>
128 double initial_time,
int size_i)
132 step =
new int[size] ;
133 step[0] = initial_step ;
134 for (
int j=1; j<size; j++) {
135 step[j] = UNDEF_STEP ;
138 the_time =
new double[size] ;
139 the_time[0] = initial_time ;
140 for (
int j=1; j<size; j++) {
141 the_time[j] = -1e20 ;
144 val =
new TyT*[size] ;
145 val[0] =
new TyT(initial_value) ;
146 for (
int j=1; j<size; j++) {
153 template<
typename TyT>
158 step =
new int[size] ;
159 for (
int j=0; j<size; j++) {
160 step[j] = UNDEF_STEP ;
163 the_time =
new double[size] ;
164 for (
int j=0; j<size; j++) {
165 the_time[j] = -1e20 ;
168 val =
new TyT*[size] ;
169 for (
int j=0; j<size; j++) {
177 template<
typename TyT>
180 pos_jtop(evo.pos_jtop) {
182 step =
new int[size] ;
183 for (
int j=0; j<size; j++) {
184 step[j] = evo.
step[j] ;
187 the_time =
new double[size] ;
188 for (
int j=0; j<size; j++) {
192 val =
new TyT*[size] ;
193 for (
int j=0; j<size; j++) {
194 if (evo.
val[j] != 0x0) {
195 val[j] =
new TyT( *(evo.
val[j]) ) ;
213 template<
typename TyT>
219 for (
int j=0; j<size; j++) {
220 if (val[j] != 0x0)
delete val[j] ;
234 template<
typename TyT>
237 cerr <<
"void Evolution<TyT>::operator= : not implemented yet ! \n" ;
243 template<
typename TyT>
246 if ( !(is_known(j)) ) return ;
249 int pos = position(j) ;
251 assert( val[pos] != 0x0) ;
255 step[pos] = UNDEF_STEP ;
256 the_time[pos] = -1e20 ;
258 if (pos == pos_jtop) {
260 while ( (pos_jtop>=0) && (val[pos_jtop] == 0x0) ) pos_jtop-- ;
272 template<
typename TyT>
275 assert(pos_jtop >= 0) ;
276 int jmax = step[pos_jtop] ;
278 if (j == jmax)
return pos_jtop ;
282 if ( (j>=step[0]) && (j<jmax) ) {
284 for (
int i=pos_jtop-1; i>=0; i--) {
293 cerr <<
"Evolution<TyT>::position: time step j = " <<
294 j <<
" not found !" << endl ;
301 template<
typename TyT>
304 assert( (j>=0) && (j<=pos_jtop) ) ;
305 assert( step[j] != UNDEF_STEP ) ;
309 while ( (n_pos == -1) && ( j < pos_jtop ) ) {
312 if (step[j] != UNDEF_STEP) n_pos = j ;
317 template<
typename TyT>
320 assert( (j>=0) && (j<=pos_jtop) ) ;
321 assert( step[j] != UNDEF_STEP ) ;
325 while ( (n_pos == -1) && ( j > 0 ) ) {
328 if (step[j] != UNDEF_STEP) n_pos = j ;
334 template<
typename TyT>
337 if (pos_jtop == -1)
return false ;
339 assert(pos_jtop >= 0) ;
341 int jmax = step[pos_jtop] ;
344 return ( val[pos_jtop] != 0x0 ) ;
347 if ( (j>=step[0]) && (j<jmax) ) {
349 for (
int i=pos_jtop-1; i>=0; i--) {
351 if (step[i] == j)
return ( val[i] != 0x0 ) ;
359 template<
typename TyT>
362 TyT* pval = val[position(j)] ;
363 assert(pval != 0x0) ;
369 template<
typename TyT>
372 int imin = position( j_min() ) ;
374 if ( ( t < the_time[ imin ] ) || ( t > the_time[pos_jtop] ) ) {
375 cerr <<
"Evolution<TyT>::operator()(double t, int order) : \n" 376 <<
"Requested time outside stored range!" << endl ;
379 assert ( order <= 2 ) ;
380 assert ( pos_jtop > order ) ;
384 while ((the_time[j0] < t) && ( j0 < pos_jtop )) {
385 j0 = next_position(j0) ;
400 int j1 = ( (j0 > imin) ? previous_position(j0) : next_position(j0) ) ;
403 double x0 = the_time[j1] - t ;
404 double x1 = the_time[j0] - t ;
405 double dx = the_time[j0] - the_time[j1] ;
407 double a0 = x1 / dx ;
408 double a1 = x0 / dx ;
410 return (a0*(*val[j0]) - a1*(*val[j1])) ;
417 int j1 = ( (j0 == imin) ? next_position(j0) : previous_position(j0) ) ;
422 j2 = next_position(j1) ;
425 j2 = next_position(j0) ;
427 j2 = previous_position(j1) ;
431 double x0 = t - the_time[j0] ;
432 double x1 = t - the_time[j1] ;
433 double dx = the_time[j0] - the_time[j1] ;
435 double x2 = t - the_time[j2] ;
437 double dx0 = the_time[j1] - the_time[j2] ;
438 double dx1 = the_time[j0] - the_time[j2] ;
441 double a0 = ( x2*x1 ) / ( dx2*dx1 ) ;
442 double a1 = ( x0*x2 ) / ( dx0*dx2 ) ;
443 double a2 = ( x0*x1 ) / ( dx0*dx1 ) ;
445 return ( a0*(*val[j0]) - a1*(*val[j1]) + a2*(*val[j2]) ) ;
451 cerr <<
" Evolution<TyT>::operator()(double t, int order) : \n" << endl ;
452 cerr <<
" The case order = " << order <<
" is not implemented!" << endl ;
462 template<
typename TyT>
465 int resu = UNDEF_STEP ;
466 for (
int i=0; i<=pos_jtop; i++) {
467 if (step[i] != UNDEF_STEP) {
473 if (resu == UNDEF_STEP) {
474 cerr <<
"Evolution<TyT>::j_min() : no valid time step found !" << endl ;
481 template<
typename TyT>
484 if (pos_jtop == -1) {
485 cerr <<
"Evolution<TyT>::j_max() : no valid time step found !" << endl ;
489 assert(pos_jtop >=0) ;
490 int jmax = step[pos_jtop] ;
491 assert(jmax != UNDEF_STEP) ;
503 template<
typename TyT>
507 TyT resu (
operator[](j) ) ;
516 int pos = position(j) ;
518 assert ( step[pos-1] != UNDEF_STEP ) ;
524 double dt = the_time[pos] - the_time[pos-1] ;
526 return ( (*val[pos]) - (*val[pos-1]) ) / dt ;
533 assert ( step[pos-2] != UNDEF_STEP ) ;
534 double dt = the_time[pos] - the_time[pos-1] ;
536 double dt2 = the_time[pos-1] - the_time[pos-2] ;
537 if (fabs(dt2 -dt) > 1.e-13) {
539 "Evolution<TyT>::time_derive: the current version is" 540 <<
" valid only for \n" 541 <<
" a constant time step !" << endl ;
545 return ( 1.5* (*val[pos]) - 2.* (*val[pos-1]) + 0.5* (*val[pos-2]) ) / dt ;
552 assert ( step[pos-2] != UNDEF_STEP ) ;
553 assert ( step[pos-3] != UNDEF_STEP ) ;
554 double dt = the_time[pos] - the_time[pos-1] ;
556 double dt2 = the_time[pos-1] - the_time[pos-2] ;
557 double dt3 = the_time[pos-2] - the_time[pos-3] ;
558 if ((fabs(dt2 -dt) > 1.e-13)||(fabs(dt3 -dt2) > 1.e-13)) {
560 "Evolution<TyT>::time_derive: the current version is valid only for \n" 561 <<
" a constant time step !" << endl ;
565 return ( 11.*(*val[pos]) - 18.*(*val[pos-1]) + 9.*(*val[pos-2])
566 - 2.*(*val[pos-3])) / (6.*dt) ;
570 cerr <<
"Evolution<TyT>::time_derive: the case n = " << n
571 <<
" is not implemented !" << endl ;
578 return operator[](j) ;
588 template<
typename TyT>
591 ofstream fich(filename) ;
593 time_t temps = time(0x0) ;
595 fich <<
"# " << filename <<
" " << ctime(&temps) ;
596 fich <<
"# " << size <<
" size" <<
'\n' ;
597 fich <<
"# " << pos_jtop <<
" pos_jtop" <<
'\n' ;
598 fich <<
"# t value... \n" ;
601 fich.setf(ios::scientific) ;
603 for (
int i=0; i<=pos_jtop; i++) {
604 if (step[i] != UNDEF_STEP) {
605 fich << the_time[i] ; fich.width(23) ;
606 assert(val[i] != 0x0) ;
607 write_formatted(*(val[i]), fich) ;
TyT ** val
Array of pointers onto the values (size = size).
Evolution(const TyT &initial_value, int initial_j, double initial_time, int initial_size)
Constructor from some initial value.
int * step
Array of time step indices (size = size).
double * the_time
Array of values of t at the various time steps (size = size).
Time evolution (*** under development ***).