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
|
+
+
+
|
#! /bin/sh
######################################
##### Model Transformation Tools #####
######################################
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.62 2001/08/01 22:14:32 geraint
## Bug fix for dassl.
##
## Revision 1.61 2001/08/01 04:06:07 geraint
## Added -i dassl for -cc and -oct.
##
## Revision 1.60 2001/07/16 22:23:00 geraint
## Fixed misleading variable name in .cc rep.
##
## Revision 1.59 2001/07/13 04:54:04 geraint
|
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
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
|
588
589
590
591
592
593
594
595
596
597
598
599
600
601
|
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
|
#else
static octave_value_list args, f;
f = feval ("${sys}_numpar", args, 1);
return f(0).${vector_value} ();
#endif
}
#ifdef STANDALONE
ColumnVector
Fmtt_residual (const ColumnVector &X,
const ColumnVector &DX,
double t)
{
#else // !STANDALONE
DEFUN_DLD (mtt_residual, args, , "mtt_residual")
{
static ColumnVector X (MTTNX+MTTNYZ);
static ColumnVector DX (MTTNX+MTTNYZ);
static double t;
X = args(0).${vector_value} ();
DX = args(1).${vector_value} ();
t = args(2).double_value ();
#endif // STANDALONE
static ColumnVector residual (MTTNX+MTTNYZ);
static ColumnVector U (MTTNU+MTTNYZ);
static ColumnVector u (MTTNU);
static ColumnVector y (MTTNY,0.0);
static ColumnVector par (MTTNPAR);
static ColumnVector dx(MTTNX);
static ColumnVector yz(MTTNYZ);
static ColumnVector x (MTTNX);
static ColumnVector ui (MTTNYZ);
static octave_value_list new_args;
x = X.extract (0,MTTNX-1);
if (MTTNYZ > 0)
ui = X.extract (MTTNX,MTTNX+MTTNYZ-1);
#ifdef STANDALONE
par = F${sys}_numpar();
u = F${sys}_input(x,y,t,par);
#else
par = feval ("${sys}_numpar", new_args, 1)(0).${vector_value} ();
new_args(0) = octave_value (x);
new_args(1) = octave_value (u);
new_args(2) = octave_value (t);
new_args(3) = octave_value (par);
u = feval ("${sys}_input", new_args, 1)(0).${vector_value} ();
#endif
U.insert (u,0);
if (MTTNYZ > 0)
U.insert (ui,MTTNU);
#ifdef STANDALONE
dx = F${sys}_${ode} (x,U,t,par);
yz = F${sys}_ae (x,U,t,par);
#else
new_args(1) = octave_value (U);
dx = feval ("${sys}_${ode}", new_args, 1)(0).${vector_value} ();
yz = feval ("${sys}_ae", new_args, 1)(0).${vector_value} ();
#endif
for (register int i = 0; i < MTTNX; i++)
residual (i) = dx(i) - DX(i);
if (MTTNYZ > 0)
residual.insert (yz,MTTNX);
#ifdef STANDALONE
return residual;
#else // !STANDALONE
return octave_value (residual);
#endif // STANDALONE
}
inline Octave_map
mtt_simpar (void)
{
#ifdef STANDALONE
return F${sys}_simpar ();
#else
static octave_value_list args;
|
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
|
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
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
|
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
|
args (6) = octave_value (static_cast<double> (nx));
args (7) = octave_value (static_cast<double> (nyz));
args (8) = octave_value (open_switches);
f = feval ("mtt_dassl", args, 1);
return f(0).${vector_value} ();
#endif
}
#ifdef STANDALONE
ColumnVector
Fmtt_residual (const ColumnVector &X,
const ColumnVector &DX,
double t)
{
#else // !STANDALONE
DEFUN_DLD (mtt_residual, args, , "mtt_residual")
{
static ColumnVector X (MTTNX+MTTNYZ);
static ColumnVector DX (MTTNX+MTTNYZ);
static double t;
X = args(0).${vector_value} ();
DX = args(1).${vector_value} ();
t = args(2).double_value ();
#endif // STANDALONE
static ColumnVector residual (MTTNX+MTTNYZ);
static ColumnVector U (MTTNU+MTTNYZ);
static ColumnVector u (MTTNU);
static ColumnVector y (MTTNY,0.0);
static ColumnVector par (MTTNPAR);
static ColumnVector dx(MTTNX);
static ColumnVector yz(MTTNYZ);
static ColumnVector x (MTTNX);
static ColumnVector ui (MTTNYZ);
static octave_value_list new_args;
x = X.extract (0,MTTNX-1);
if (MTTNYZ > 0)
ui = X.extract (MTTNX,MTTNX+MTTNYZ-1);
#ifdef STANDALONE
par = F${sys}_numpar();
u = F${sys}_input(x,y,t,par);
#else
par = feval ("${sys}_numpar", new_args, 1)(0).${vector_value} ();
new_args(0) = octave_value (x);
new_args(1) = octave_value (u);
new_args(2) = octave_value (t);
new_args(3) = octave_value (par);
u = feval ("${sys}_input", new_args, 1)(0).${vector_value} ();
#endif
U.insert (u,0);
if (MTTNYZ > 0)
U.insert (ui,MTTNU);
#ifdef STANDALONE
dx = F${sys}_${ode} (x,U,t,par);
yz = F${sys}_ae (x,U,t,par);
#else
new_args(1) = octave_value (U);
dx = feval ("${sys}_${ode}", new_args, 1)(0).${vector_value} ();
yz = feval ("${sys}_ae", new_args, 1)(0).${vector_value} ();
#endif
for (register int i = 0; i < MTTNX; i++)
residual (i) = dx(i) - DX(i);
if (MTTNYZ > 0)
residual.insert (yz,MTTNX);
#ifdef STANDALONE
return residual;
#else // !STANDALONE
return octave_value (residual);
#endif // STANDALONE
}
EOF
;;
"euler" | "rk4" | *)
cat <<EOF >> $filename
inline ColumnVector
mtt_euler (ColumnVector &x,
const ColumnVector &dx,
|