Annotate Patch Side-by-side Diff

Differences From Artifact [a0ec3f6c8d]:

To Artifact [edf3325ddc]:

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
31
32
33
34
35
36












+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+








#! /bin/sh

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

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.57.2.4  2001/07/02 00:34:56  geraint
## gcc-3.0 compatibility.
##
## Revision 1.57.2.3  2001/06/25 23:28:29  geraint
## Generic mtt_rate and mtt_output - allows method independent calls.
##
## 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

235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254








+
+
+
+
+
+








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

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
395
396
397
398
399
400








-
-
-
+
+
+









+
+



+
+
+
+
+
+








	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"

#include "mtt_${algebraic_solver}.hh"

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

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
451

452








453
454
455
456
457
458
459

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
493
494
495
496
497
498

499
500
501
502
503
504
505
506

507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523








+
+
+


+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+







+
+
+

-
+







-
+

+
+
+
+
+
+
+
+








	ColumnVector &u,
	const double &t,
	ColumnVector &par);
 
EOF
fi
cat <<EOF >> $filename

void set_signal_handlers (void);

#endif // STANDALONE

ColumnVector
mtt_ae (ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par)
{
#ifdef STANDALONE
    return F${sys}_ae(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}_ae", args, 1);
    return f(0).${vector_value} ();
#endif
}

inline ColumnVector
mtt_input (ColumnVector &x,
	   ColumnVector &y,
	   const double &t,
	   ColumnVector &par)
{
  static ${algebraic_solver} ae(mtt_ae,MTTNPAR,MTTNU,MTTNX,MTTNY,MTTNYZ);
  static ColumnVector u  (MTTNU);

#ifdef STANDALONE
  return F${sys}_input (x, y, t, par);
  u = 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);
  return f(0).${vector_value} ();
  u = f(0).${vector_value} ();
#endif
  if (MTTNYZ == 0)
    {
      return u;
    }
  else
    {
      return ae.solve(x,u,t,par);
    }
}

inline ColumnVector
mtt_logic (ColumnVector &x,
	   ColumnVector &u,
	   const double &t,
	   ColumnVector &par)
︙
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

577
578
579
580
581
582
583




584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602

603
604
605
606
607
608
609
610








-
-
-
-
+
+
+
+















-
+








  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
︙
650
651
652
653
654
655
656
657

658
659
660
661
662
663


664
665



666
667
668
669
670
671
672

714
715
716
717
718
719
720

721
722
723
724
725
726
727
728
729


730
731
732
733
734
735
736
737
738
739








-
+






+
+
-
-
+
+
+









inline void
mtt_write (const double &t,
	   ColumnVector &x,
	   ColumnVector &y,
	   const int &nrows,
	   const bool dump_data = false,
	   ostream &file = cout)
	   std::ostream &file = std::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)
︙
680
681
682
683
684
685
686
687
688


689
690
691
692
693
694
695
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
723
724

725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744

745
746
747
748
749
750
751
752
753

747
748
749
750
751
752
753


754
755
756
757
758
759
760
761
762
763

764
765

766
767
768
769


770
771
772


773
774
775
776
777
778
779

780
781
782
783
784
785
786
787
788

789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808

809


810
811
812
813
814
815
816








-
-
+
+








-
+

-
+



-
-



-
-
+
+





-
+








-
+



















-
+
-
-








    data.elem (row, ++col) = x.elem (i);

  row++;

  if (nrows == row)
    {
#ifdef STANDALONE
      $save_ascii_data_function (file, data, "MTT_data");
//    cout << data << endl;
      $save_ascii_data_function (file, data, "mtt_data");
//    std::cout << data << std::endl;
#else // ! STANDALONE
  set_global_value ("MTT_data", data);
#endif
    row = 0;
    }
}

#ifdef STANDALONE
void dump_data (ostream &file)
void dump_data (std::ostream &file)
{
  ColumnVector null (0.0);
  ColumnVector null (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");
  std::cerr << "# Writing data to MTT.core (signal " << signum << ")" << std::endl;
  std::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;
      default:
        cerr << "# Warning: make_ode2odes needs updating!" << endl;
        std::cerr << "# Warning: make_ode2odes needs updating!" << std::endl;
	signal (signum, SIG_DFL);
	raise (signum);
        break;
    }
  corefile.close ();
  set_signal_handlers ();
}

void set_signal_handlers (void)
{
  signal (SIGFPE,	handle_signal);
  signal (SIGINT,	handle_signal);
  signal (SIGQUIT,	handle_signal);
}

int main (void) {
  set_signal_handlers ();
#else
DEFUN_DLD (${sys}_ode2odes, args, ,
"Octave ode2odes representation of system with $method integration method
"Octave ode2odes representation of system with $method integration method\nUsage: mtt_data = ${sys}_ode2odes (x0, par, simpar)\n")
Usage: mtt_data = ${sys}_ode2odes (x0, par, simpar)
")
{
  static octave_value_list retval;
#endif // STANDALONE
  static ColumnVector	x0;
  static ColumnVector	par;
  static Octave_map	simpar;
︙
824
825
826
827
828
829
830
831

832
833
834
835
836
837
838

887
888
889
890
891
892
893

894
895
896
897
898
899
900
901








-
+








    {
      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

911
912
913
914
915
916
917




918
919
920
921
922
923
924
925
926

927
928
929
930
931
932
933
934








-
-
-
-
+
+
+
+





-
+








          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 ]