Annotate Patch Side-by-side Diff

Differences From Artifact [a0ec3f6c8d]:

To Artifact [d8ff4bb840]:

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












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








#! /bin/sh

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

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## 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

232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251








+
+
+
+
+
+








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

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
395
396
397








-
-
-
+
+
+












+
+
+
+
+
+
+
+








	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

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
493
494
495








-
+
+

+








+
+
+
+
+
+
-
+
+
+
+
+
+
+
+
+








	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

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
589
590
591








-
-
-
-
+
+
+
+















-
+








  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

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








-
+






+
+
-
-
+
+
+









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

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


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








-
+








-
+

-
+



-
-



-
-
+
+





-
+








-
+



















-
+
-
-









  row++;

  if (nrows == row)
    {
#ifdef STANDALONE
      $save_ascii_data_function (file, data, "MTT_data");
//    cout << data << endl;
//    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

868
869
870
871
872
873
874

875
876
877
878
879
880
881
882








-
+








    {
      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

892
893
894
895
896
897
898




899
900
901
902
903
904
905
906
907

908
909
910
911
912
913
914
915








-
-
-
-
+
+
+
+





-
+








          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 ]