Differences From Artifact [a0ec3f6c8d]:

To Artifact [1f5fcd3601]:


1
2
3
4
5
6
7
8
9
10
11












12
13
14
15
16
17
18
#! /bin/sh

     ###################################### 
     ##### Model Transformation Tools #####
    ######################################

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$












## Revision 1.56  2001/03/30 15:13:58  gawthrop
## Rationalised simulation modes to each return mtt_data
##
## Revision 1.55  2001/03/27 13:21:59  geraint
## Octave version compatibility for save_ascii_data(_for_plotting).
##
## Revision 1.54  2001/03/27 01:14:27  geraint











>
>
>
>
>
>
>
>
>
>
>
>







1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#! /bin/sh

     ###################################### 
     ##### Model Transformation Tools #####
    ######################################

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.57.2.2  2001/06/05 03:20:40  geraint
## added -ae option to select algebraic equation solution method.
##
## Revision 1.57.2.1  2001/05/04 04:07:24  geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.57  2001/04/01 03:38:54  geraint
## Reset row to zero after write to file, ready for subsequent runs.
## Eliminates SIGSEGV in Octave when _ode2odes called multiple times.
##
## Revision 1.56  2001/03/30 15:13:58  gawthrop
## Rationalised simulation modes to each return mtt_data
##
## Revision 1.55  2001/03/27 13:21:59  geraint
## Octave version compatibility for save_ascii_data(_for_plotting).
##
## Revision 1.54  2001/03/27 01:14:27  geraint
217
218
219
220
221
222
223






224
225
226
227
228
229
230
filename=${sys}_ode2odes.${lang}

if [ -n "$3" ]; then
  method=$3    
else
  method=implicit  
fi







echo Creating $filename with $method integration method

# Find system constants
Nx=`mtt_getsize $sys x` # States
Nu=`mtt_getsize $sys u` # Inputs 
Ny=`mtt_getsize $sys y` # Inputs  







>
>
>
>
>
>







229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
filename=${sys}_ode2odes.${lang}

if [ -n "$3" ]; then
  method=$3    
else
  method=implicit  
fi

if [ -n "$4" ]; then
    algebraic_solver=$4
else
  algebraic_solver="Reduce_Solver"
fi

echo Creating $filename with $method integration method

# Find system constants
Nx=`mtt_getsize $sys x` # States
Nu=`mtt_getsize $sys u` # Inputs 
Ny=`mtt_getsize $sys y` # Inputs  
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361








362
363
364
365
366
367
368
	feval_header=parse.h
	save_ascii_data_function=save_ascii_data_for_plotting
	;;
esac

cat <<EOF  > $filename
#include <octave/oct.h>
#include <octave/ov-struct.h>
#include <octave/load-save.h>
#include <octave/lo-mappers.h>
#include <octave/variables.h>

#ifndef STANDALONE
#include <octave/${feval_header}>
#endif

#include "${sys}_def.h"
#include "${sys}_sympar.h"

#ifdef STANDALONE
#include <csignal>
#include <fstream>









extern ColumnVector F${sys}_input (
	ColumnVector &x,
	ColumnVector &y,
	const double &t,
	ColumnVector &par);








|
|
|












>
>
>
>
>
>
>
>







358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
	feval_header=parse.h
	save_ascii_data_function=save_ascii_data_for_plotting
	;;
esac

cat <<EOF  > $filename
#include <octave/oct.h>
#include <octave/load-save.h>
#include <octave/lo-mappers.h>
#include <octave/ov-struct.h>
#include <octave/variables.h>

#ifndef STANDALONE
#include <octave/${feval_header}>
#endif

#include "${sys}_def.h"
#include "${sys}_sympar.h"

#ifdef STANDALONE
#include <csignal>
#include <fstream>

#include "mtt_${algebraic_solver}.hh"

extern ColumnVector F${sys}_ae (
	ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par);

extern ColumnVector F${sys}_input (
	ColumnVector &x,
	ColumnVector &y,
	const double &t,
	ColumnVector &par);

426
427
428
429
430
431
432
433

434

435
436
437
438
439
440
441
442






443








444
445
446
447
448
449
450
	ColumnVector &u,
	const double &t,
	ColumnVector &par);
 
EOF
fi
cat <<EOF >> $filename
#endif // STANDALONE




inline ColumnVector
mtt_input (ColumnVector &x,
	   ColumnVector &y,
	   const double &t,
	   ColumnVector &par)
{
#ifdef STANDALONE






  return F${sys}_input (x, y, t, par);








#else
  static octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (y);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
  f = feval ("${sys}_input", args, 1);







|
>

>








>
>
>
>
>
>
|
>
>
>
>
>
>
>
>







452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
	ColumnVector &u,
	const double &t,
	ColumnVector &par);
 
EOF
fi
cat <<EOF >> $filename

void set_signal_handlers (void);

#endif // STANDALONE

inline ColumnVector
mtt_input (ColumnVector &x,
	   ColumnVector &y,
	   const double &t,
	   ColumnVector &par)
{
#ifdef STANDALONE
  static ColumnVector u  (MTTNU);
  static ColumnVector ui (MTTNYZ);
  static ColumnVector U  (MTTNU+MTTNYZ);
 
  static ${algebraic_solver} ae(F${sys}_ae,MTTNPAR,MTTNU,MTTNX,MTTNY,MTTNYZ);

  u = F${sys}_input (x, y, t, par);
  if (MTTNYZ == 0)
    {
      return u;
    }
  else
    {
      return ae.solve(x,u,t,par);
    }
#else
  static octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (y);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
  f = feval ("${sys}_input", args, 1);
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
  args (0) = octave_value (x);
  f = feval ("${sys}_state", args, 1);
  return f(0).${vector_value} ();
#endif
}

inline ColumnVector
mtt_${ode} (ColumnVector &x,
	    ColumnVector &u,
	    const double &t,
	    ColumnVector &par)
{
#ifdef STANDALONE
  return F${sys}_${ode} (x, u, t, par);
#else
  static octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (u);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
  f = feval ("${sys}_${ode}", args, 1);
  return f(0).${vector_value} ();
#endif
}

inline ColumnVector
mtt_${odeo} (ColumnVector &x,
	     ColumnVector &u,
	     const double &t,
	     ColumnVector &par)
{
#ifdef STANDALONE
  return F${sys}_${odeo} (x, u, t, par);
#else







|
|
|
|















|







555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
  args (0) = octave_value (x);
  f = feval ("${sys}_state", args, 1);
  return f(0).${vector_value} ();
#endif
}

inline ColumnVector
mtt_rate (ColumnVector &x,
	  ColumnVector &u,
	  const double &t,
	  ColumnVector &par)
{
#ifdef STANDALONE
  return F${sys}_${ode} (x, u, t, par);
#else
  static octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (u);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
  f = feval ("${sys}_${ode}", args, 1);
  return f(0).${vector_value} ();
#endif
}

inline ColumnVector
mtt_output (ColumnVector &x,
	     ColumnVector &u,
	     const double &t,
	     ColumnVector &par)
{
#ifdef STANDALONE
  return F${sys}_${odeo} (x, u, t, par);
#else
657
658
659
660
661
662
663


664
665

666
667
668
669
670
671
672
	   ostream &file = cout)
{
  static Matrix data;
  static int row;

  if (dump_data)
    {


       Matrix written_data = data.extract (0, 0, row-1, data.cols ()-1);
       $save_ascii_data_function (file, written_data, "MTT_data");

       return;
    }

  const int nx = x.length (), ny = y.length ();
  register int col = 0;

  if (0 == row)







>
>
|
|
>







699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
	   ostream &file = cout)
{
  static Matrix data;
  static int row;

  if (dump_data)
    {
       if (row > 0)
         {
           Matrix written_data = data.extract (0, 0, row-1, data.cols ()-1);
           $save_ascii_data_function (file, written_data, "MTT_data");
	 }
       return;
    }

  const int nx = x.length (), ny = y.length ();
  register int col = 0;

  if (0 == row)
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
#ifdef STANDALONE
void dump_data (ostream &file)
{
  ColumnVector null (0.0);
  mtt_write (0.0, null, null, 0, true, file);
}

void set_signal_handlers (void);

void handle_signal (int signum)
{
  // handle some signals to ensure data is written.
  cerr << "# Writing data to MTT.core (signal " << signum << ")" << endl;
  ofstream corefile ("MTT.core");
  dump_data (corefile);
  switch (signum)
    {
      case SIGFPE:
        // Intel chips do not raise SIGFPE for DIVZERO :-(
        raise (SIGABRT);
        break;
      case SIGINT:
	break;
      case SIGQUIT:
	signal (SIGQUIT, SIG_DFL);
	raise (SIGQUIT);
	break;







<
<










|







741
742
743
744
745
746
747


748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
#ifdef STANDALONE
void dump_data (ostream &file)
{
  ColumnVector null (0.0);
  mtt_write (0.0, null, null, 0, true, file);
}



void handle_signal (int signum)
{
  // handle some signals to ensure data is written.
  cerr << "# Writing data to MTT.core (signal " << signum << ")" << endl;
  ofstream corefile ("MTT.core");
  dump_data (corefile);
  switch (signum)
    {
      case SIGFPE:
        // Intel chips do not raise SIGFPE for DIVZERO :-(
//        raise (SIGABRT);
        break;
      case SIGINT:
	break;
      case SIGQUIT:
	signal (SIGQUIT, SIG_DFL);
	raise (SIGQUIT);
	break;
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
    {
      x (i) = x0 (i);
    }

  for (register int j = 0, i = 1; i <= ilast; i++)
    {
      u	= mtt_input (x, y, t, par);
      y	= mtt_${odeo} (x, u, t, par);
      if (0 == j)
	{
           mtt_write (t, x, y, nrows);
	}
EOF
if [ "$method" = "rk4" ]; then
cat << EOF >> $filename







|







867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
    {
      x (i) = x0 (i);
    }

  for (register int j = 0, i = 1; i <= ilast; i++)
    {
      u	= mtt_input (x, y, t, par);
      y	= mtt_output (x, u, t, par);
      if (0 == j)
	{
           mtt_write (t, x, y, nrows);
	}
EOF
if [ "$method" = "rk4" ]; then
cat << EOF >> $filename
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
          t2 = t + ddt;

        ColumnVector
          x1 (x),
          x2 (x),
          x3 (x);

        k1 = ddt * mtt_${ode} (x , u, t , par); x1 += k1 * 0.5;
        k2 = ddt * mtt_${ode} (x1, u, t1, par); x2 += k2 * 0.5;
        k3 = ddt * mtt_${ode} (x2, u, t1, par); x3 += k3;
        k4 = ddt * mtt_${ode} (x3, u, t2, par);
        dx = (k1 + 2.0 * (k2 + k3) + k4) / (6.0 * ddt);
      }
EOF
else
cat << EOF >> $filename
      dx = mtt_${ode} (x, u, t, par);
EOF
fi

if [ "$method" = "implicit" ]; then
cat <<EOF >> $filename

      AA = mtt_smxa (x, u, ddt, par);







|
|
|
|





|







891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
          t2 = t + ddt;

        ColumnVector
          x1 (x),
          x2 (x),
          x3 (x);

        k1 = ddt * mtt_rate (x , u, t , par); x1 += k1 * 0.5;
        k2 = ddt * mtt_rate (x1, u, t1, par); x2 += k2 * 0.5;
        k3 = ddt * mtt_rate (x2, u, t1, par); x3 += k3;
        k4 = ddt * mtt_rate (x3, u, t2, par);
        dx = (k1 + 2.0 * (k2 + k3) + k4) / (6.0 * ddt);
      }
EOF
else
cat << EOF >> $filename
      dx = mtt_rate (x, u, t, par);
EOF
fi

if [ "$method" = "implicit" ]; then
cat <<EOF >> $filename

      AA = mtt_smxa (x, u, ddt, par);

MTT: Model Transformation Tools
GitHub | SourceHut | Sourceforge | Fossil RSS ]