Annotate Patch Side-by-side Diff

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

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

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

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/ov-struct.h>
#include <octave/load-save.h>
#include <octave/lo-mappers.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

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
#endif // STANDALONE

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);

  return F${sys}_input (x, y, t, par);
  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

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_${ode} (ColumnVector &x,
	    ColumnVector &u,
	    const double &t,
	    ColumnVector &par)
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_${odeo} (ColumnVector &x,
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

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");
           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

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 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);
//        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

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_${odeo} (x, u, 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

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_${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);
        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_${ode} (x, u, t, par);
      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 ]