Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-05 - 43,50337/21/cadsim.sim
There is 1 other file named cadsim.sim in the archive. Click here to see a list.
00100	OPTIONS(/E);
00200	EXTERNAL PROCEDURE abort;
00300	EXTERNAL REF(Infile) PROCEDURE findinfile;
00400	EXTERNAL REF(Outfile) PROCEDURE findoutfile;
00500	EXTERNAL TEXT PROCEDURE conc,upcase,frontstrip,rest,
00600	checkextension;
00700	EXTERNAL CHARACTER PROCEDURE fetchar,findtrigger;
00800	EXTERNAL LONG REAL PROCEDURE scanreal;
00900	EXTERNAL INTEGER PROCEDURE checkreal,checkint,scanint,ilog;
01000	EXTERNAL BOOLEAN PROCEDURE menu;
01100	EXTERNAL CLASS simeio;
01200	COMMENT----------------------------------------------------------------+
01300	!                                                                      !
01400	!                                CADSIM                                !
01500	!                                ------                                !
01600	!                                                                      !
01700	!                 CONTINUOUS AND DISCRETE SIMULATION                   !
01800	!                                                                      !
01900	!         BY ROBERT J. W. SIM                                          !
02000	!            DEPARTMENT OF COMPUTING AND CONTROL,                      !
02100	!            IMPERIAL COLLEGE,                                         !
02200	!            LONDON S.W.7                                              !
02300	!                                                                      !
02400	!            IMPLEMENTED JULY 1975                                     !
02500	!                                                                      !
02600	!   THIS VERSION NO 1.0 1977-12-13 BY                                  !
02700	!            G.LYMAN AND I.NORDSTROEM                                  !
02800	!            DEPARTMENT OF WEAPON SYSTEMS                              !
02900	!                                                                      !
03000	!            M.OHLIN                                                   !
03100	!            DEPARTMENT OF PLANNING AND OPERATIONS RESEARCH            !
03200	!                                                                      !
03300	!            SWEDISH NATIONAL DEFENCE RESEARCH INSTITUTE               !
03400	!            S-104 50 STOCKHOLM 80                                     !
03500	!      							               !
03600	!   CADSIM IS A SUBCLASS OF SIMULATION AND SAFEIO DESIGNED TO PROVIDE  !
03700	!   FACILITIES FOR COMBINED DISCRETE EVENT AND CONTINUOUS SYSTEM       !
03800	!   SIMULATION.                                                        !
03900	!   CADSIM HAS ONE CLASS ATTRIBUTE                                     !
04000	!                                                                      !
04100	!           1) CLASS CONTINUOUS --- THIS IS THE COUNTERPART OF CLASS   !
04200	!                                  PROCESS IN DISCRETE EVENT SIMULATION!
04300	!                                  IT IS A SUBCLASS OF PROCESS.        !
04400	!                                                                      !
04500	+----------------------------------------------------------------------+
04600	;
04700	simeio CLASS cadsim;
04800	
04900	BEGIN
05000	
05100	    OPTIONS(/P:" ---*** CONTINUOUS ***---");
05200	
05300	    COMMENT----------------------------------------------------------------+
05400	    !                                                                      !
05500	    !   CLASS CONTINUOUS PROVIDES THE FRAMEWORK FOR OBJECTS WHICH OBEY     !
05600	    !   DIFFERENTIAL DYNAMICS.                                             !
05700	    !                                                                      !
05800	    !                                                                      !
05900	    !   THE PARAMETER, N, IS USED TO SPECIFY THE NUMBER OF FIRST ORDER     !
06000	    !   EQUATIONS USED IN PROCEDURE DERIVATIVES. ITS VALUE IS USED TO      !
06100	    !   DIMENSION ALL USER ACCESSIBLE ARRAYS AND IN PROCEDURES ERROR       !
06200	    !   AND TAKE_A_STEP.                                                   !
06300	    !                                                                      !
06400	    !                         VIRTUAL PROCEDURES                           !
06500	    !                                                                      !
06600	    !   DERIVATIVES   :THIS PROCEDURE GIVES THE RELATIONSHIP BETWEEN THE   !
06700	    !                  ELEMENTS OF THE RATE AND SATE VECTORS.              !
06800	    !                  IT MUST ALWAYS BE WRITTEN BY THE USER.              !
06900	    !                                                                      !
07000	    !   ERROR         :THE USER MAY REWRITE THIS IF HE BELIEVES THAT HE    !
07100	    !                  NEEDS A DIFFERENT WAY OF OBTAINING THE CRITERIA FOR !
07200	    !                  HALVING THE STEP LENGTH.                            !
07300	    !                                                                      !
07400	    !   TAKE_A_STEP   :THIS MAY BE REWRITTEN IF THE USER REQUIRES A        !
07500	    !                  DIFFERENT METHOD OF ADVANCING THE INTEGRATION BY    !
07600	    !                  ONE STEP.                                           !
07700	    !                                                                      !
07800	    !   THE VARIABLES LISTED BELOW ARE GLOBAL TO CONTINUOUS OBJECTS AND    !
07900	    !   ARE THEREFOR AVAILABLE TO THE USER.                                !
08000	    !                                                                      !
08100	    !   MAX_H         :USER DEFINED MAXIMUM INTEGRATION STEP LENGTH        !
08200	    !   MIN_H         :USER DEFINED MINIMUM INTEGRATION STEP LENGTH        !
08300	    !   ERROR_LIMIT   :THE ACCURACY WHICH THE USER REQUIRES FROM THE       !
08400	    !                  INTEGRATION ROUTINE.                                !
08500	    !   T             :THE TIME RELATIVE TO THE START OF INTEGRATION.      !
08600	    !                                                                      !
08700	    !   RATE          :A VECTOR WHICH HOLDS THE VALUES OF THE RATE OF      !
08800	    !                  CHANGE OF THE STATE VARIABLES.                      !
08900	    !   STATE         :A VECTOR HOLDING THE VALUE OF THE STATE VARIABLES.  !
09000	    !   ERROR_VECTOR  :A VECTOR HOLDING THE ESTIMATED ERROR OF EACH STATE  !
09100	    !                  VARIABLE.                                           !
09200	    !   INITIAL_STATE :A VECTOR WHICH HOLDS THE VALUES OF THE STATE        !
09300	    !                  VARIABLES AT THE BEGINING OF EACH TIME STEP.        !
09400	    !                                                                      !
09500	    !   SLOPE_1, SLOPE_2, SLOPE_3           : VECTORS WHICH HOLD THE VALUE ! 
09600	    !                                         OF THE SLOPE OF THE FUNCTION !
09700	    !                                         AT POINTS WITHIN THE STEP.   !
09800	    !                                         THESE ARE NECESSARY FOR THE  !
09900	    !                                          INTERPOLATION-PROCEDURES    !
10000	    !   RELATIVE      :A BOOLEAN INDICATING IF A RELATIVE OR ABSOLUTE ERROR!
10100	    !                  TEST IS TO BE MADE.                                 !
10200	    !                  I.E. RELATIVE = TRUE  ==> RELATIVE ERROR TEST       !
10300	    !                       RELATIVE = FALSE ==> ABSOLUTE ERROR TEST       !
10400	    !   STOP          :A BOOLEAN USED TO STOP THE INTEGRATION WHEN THE     !
10500	    !                  NORMAL END CONDITIONS ARE STILL TRUE C.F.'INTEGRATE'!
10600	    !   TRACE_LEVEL   :AN INTEGER DETERMINING HOW MUCH INFORMATION ABOUT   !
10700	    !                  THE INTEGRATION IS TO BE PRINTED.                   !
10800	    !                                                                      !
10900	    +----------------------------------------------------------------------+
11000	;
11100	    process CLASS nondiscrete;;
11200	
11300	    nondiscrete CLASS continuous(n);
11400	    INTEGER n;
11500	
11600	    VIRTUAL:PROCEDURE derivatives,take_a_step, error;
11700	
11800	    BEGIN
11900		REAL max_h, min_h, error_limit, t, start_time;
12000		REAL ARRAY rate, state, error_vector,
12100		initial_state, slope_1, slope_2, slope_3[1:n];
12200		BOOLEAN relative, stop, trace;
12300		TEXT trace_level;REAL ARRAY spec[1000:1005];
12400		OPTIONS(/P:" --- INTEGRATE(CONTINUOUS) ---");
12500	
12600		COMMENT----------------------------------------------------------------+
12700		!                                                                      !
12800		!   INTEGRATE IS THE PROCEDURE WHICH CONTROLS THE INTEGRATION OF STATE !
12900		!   VARIABLES. IT HAS THREE LOGICAL STEPS WHICH ARE REPEATED UNTIL     !
13000		!   ONE OF THE BOOLEAN CONDITIONS BECOMES FALSE:-                      !
13100		!                                                                      !
13200		!        1) FIND AN APPROPRIATE STEP LENGTH.                           !
13300		!        2) CALL A ROUTINE TO ADVANCE THE INTEGRATION BY THE CHOSEN    !
13400		!           STEP LENGTH.                                               !
13500		!        3) ESTIMATE IF THE ERROR IS SMALL ENOUGH.                     !
13600		!                                                                      !
13700		!                             PARAMETERS                               !
13800		!                                                                      !
13900		!   NOT_STATE_EVENT    :THIS GIVES THE CONDITION OF A STATE EVENT      !
14000		!   INTEGR_TIME      :THIS GIVES THE TIME OF INTEGRATION               !
14100		!                                                                      !
14200		!        BOTH PARAMETERS GIVE THE CONDITION FOR TERMINATION OF         !
14300		!        INTEGRATION                                                   !
14400		!                                                                      !
14500		!                             PROCEDURES                               !
14600		!                                                                      !
14700		!   REDUCE        :THIS GIVES THE CRITERIA FOR REDUCING THE STEP LENGTH!
14800		!   ERROR         :OPERATES ON ERROR_VECTOR TO PRODUCE A SCALER VALUE  !
14900		!                                                                      !
15000		!                             VARIABLES                                !
15100		!                                                                      !
15200		!   H             :THE CURRENT STEP LENGTH IN USE                      !
15300		!   LARGE_H       :EITHER HOLDS THE TIME INTERVAL TO THE NEXT DISCRETE !
15400		!                  EVENT OR, IF IT IS LESS THE USER DEFINED MAX_H      !
15500		!   SAVE_H        :SAVES THE VALUE OF H WHEN IT IS GOING TO BE REDUCED !
15600		!                  TO THE VALUE OF LARGE_H.                            !
15700		!   ERROR_VALUE   :HOLDS THE SCALAR VALUE RETURNED BY ERROR            !
15800		!   I             :FOR LOOP INDEX                                      !
15900		!   NEXT          :REFERENCE TO A PROCESS OBJECT, USED TO LOOK FOR A   !
16000		!                  DISCRETE EVENT. I.E. A PROCESS NOT IN CLASS         !
16100		!                  CONTINUOUS.                                         !
16200		!                                                                      !
16300		+----------------------------------------------------------------------+
16400	;
16500		PROCEDURE integrate(not_state_event, integr_time);
16600	
16700		NAME not_state_event;
16800		BOOLEAN not_state_event;
16900		REAL integr_time;
17000	
17100		BEGIN
17200		    REAL h, large_h, error_value, save_h;
17300		    INTEGER i, j, k;
17400		    REF(process) next; INTEGER ARRAY c[0:6];
17500	
17600		    OPTIONS(/P:" REDUCE(INTEGRATE) ---");
17700	
17800		    COMMENT----------------------------------------------------------------+
17900		    !                                                                      !
18000		    !   REDUCE DETERMINES IF THE INTEGRATION STEP SIZE SHOULD BE REDUCED   !
18100		    !                                                                      !
18200		    !   REDUCE:= 'TRUE'  IF THE ERROR IS TO GREAT AND H > 2*USERS MINIMUM  !
18300		    !                   VALUE FOR H                                        !
18400		    !                                                                      !
18500		    !   REDUCE:= 'FALSE'  IF THE ERROR IS SMALL ENOUGH OR H < 2*USERS      !
18600		    !                    MINIMUM VALUE FOR H. IN THE LATTER CASE AN ERROR  !
18700		    !                    CONDITION HAS ARISEN AND THE INTEGRATION IS HALTED!
18800		    !                                                                      !
18900		    +----------------------------------------------------------------------+
19000	;
19100		    BOOLEAN PROCEDURE reduce;
19200		    BEGIN
19300			error_value:= error;
19400			IF error_value <= error_limit
19500			THEN reduce:= FALSE ELSE
19600			IF h < 2*min_h THEN
19700			abort("? CADSIM ERROR 001"
19800			"- ERROR LIMIT NOT OBTAINABLE WITH CURRENT MIN_H")
19900			ELSE reduce:= TRUE;
20000		    END --- REDUCE ---;
20100		    OPTIONS(/P:" ERROR(INTEGRATE) ---");
20200	
20300		    COMMENT----------------------------------------------------------------+
20400		    !                                                                      !
20500		    !   ERROR PRODUCES A SCALER VALUE FROM ERROR_VECTOR.IT MAY BE REWRITTEN!
20600		    !   BY THE USER. THE FORM PROVIDED SELECTS THE LARGEST ELEMENT FROM    !
20700		    !   THE VECTOR.                                                        !
20800		    !                                                                      !
20900		    +----------------------------------------------------------------------+
21000	;
21100		    REAL PROCEDURE error;
21200		    BEGIN
21300			INTEGER i;
21400			REAL max_value;
21500	
21600			FOR i:= 1 STEP 1 UNTIL n DO
21700			IF error_vector[i] > max_value
21800			THEN max_value:= error_vector[i];
21900			error:= max_value;
22000	
22100		    END --- ERROR ---;
22200		    OPTIONS(/P:" --- TRACE HEADING ---");
22300	
22400		    PROCEDURE trace_heading;
22500		    BEGIN
22600			Outtext(" SYSTEM TIME"
22700			"     STEP LENGTH"
22800			"        LARGE_H"
22900			"    ERROR VALUE"
23000			"     NEXT DISCRETE");
23100			Outimage;
23200			Setpos(18);
23300			i:=1;
23400			WHILE c[i]>0 AND c[i]<=n DO
23500			BEGIN
23600			    Outtext("STATE[");
23700			    Outint(c[i],5);
23800			    Outtext("]   "); i:=i+1;
23900			END;
24000			WHILE c[i]>=1000 AND c[i]<= 1005 DO
24100			BEGIN
24200			    Outtext("SPEC[ ");
24300			    Outint(c[i],5);
24400			    Outtext("]   ");i:=i+1;
24500			END;
24600			Outimage;
24700			Outimage;
24800		    END *** TRACE HEADING ***;
24900		    OPTIONS(/P:" --- TRACE VALUES ---;");
25000	
25100		    PROCEDURE trace_values;
25200		    BEGIN
25300			Outfix(time+h,7,13);
25400			Outfix(h,7,15);
25500			Outfix(large_h,7,15);
25600			Outfix(error_value,7,15);
25700			Outfix(IF next=/=NONE
25800			THEN next.evtime
25900			ELSE 0.0        ,7,16);
26000			Outimage; Setpos(18);
26100			i:=1;
26200			WHILE c[i]>0 AND c[i]<=n DO
26300			BEGIN
26400			    Outfix(state[c[i]],7,15);
26500			    i:=i+1;
26600			END WHILE;
26700			WHILE c[i]>=1000 AND c[i]<=1005 DO
26800			BEGIN Outfix(spec[c[i]],7,15);
26900			    i:=i+1;
27000			END WHILE;
27100			Outimage;
27200		    END *** TRACE VALUES ***;
27300		    OPTIONS(/P:"   MAIN CODE FOR INTEGRATE");
27400	
27500		    ! SET START TIME;
27600		    start_time:= time;
27700	
27800		    IF trace_level NE "" THEN trace:=TRUE ELSE trace:=FALSE;
27900	
28000	            COMMENT----------------------------------------------------+
28100	            !   TRACE_LEVEL IS A TEXT VALUE WHICH CONTAINS INFORMATION !
28200	            !   FOR THE TRACING.                                       !
28300	            !   EXAMPLE: TRACE_LEVEL:="1 2 1003"                       !
28400	            !   MEANS THAT STATE[1],STATE[2] AND SPEC[1003] WILL BE    !
28500	            !   TRACED.                                                !
28600	            +----------------------------------------------------------+
28700	;
28800		    IF trace THEN
28900		    BEGIN TEXT t1;
29000			t1:-Blanks(trace_level.Length);
29100			t1:=trace_level;
29200			i:=1;
29300			WHILE i<7 AND t1 NE "" DO
29400			BEGIN
29500			    c[i]:=t1.Getint;
29600			    t1:-t1.Sub(t1.Pos,t1.Length-t1.Pos+1);
29700			    i:=i+1;
29800			END WHILE;
29900		    END IF;
30000	
30100		    IF max_h <= 0 OR min_h <= 0 THEN
30200		    abort("? CADSIM ERROR 002 -"
30300		    " MAX_H AND MIN_H HAVE NOT BEEN INITIALISED WITH POSITIVE "
30400		    "VALUES")
30500		    ELSE stop:= FALSE;
30600	
30700		    COMMENT *** PUT DUMMY EVENT ON SQS AND EXECUTE BOGUS STEP ***;
30800	
30900		    ACTIVATE NEW process AT start_time+integr_time;
31000		    ! ?	    HOLD(1/(10**70));  hold(0);
31100	
31200		    IF trace
31300		    THEN
31400		    BEGIN
31500			trace_heading;
31600			trace_values;
31700		    END;
31800	
31900		    large_h:= h := max_h;
32000		    IF integr_time <= 0 THEN integr_time:= &37;
32100	
32200		    WHILE (not_state_event) AND (t < integr_time) AND NOT stop DO
32300		    BEGIN
32400	
32500			COMMENT *** RECORD INITIAL VALUES OF STATE VARIABLES ***;
32600	
32700			FOR i:= 1 STEP 1 UNTIL n DO initial_state[i]:= state[i];
32800	
32900			COMMENT *** LOOK FOR DISCRETE EVENT IN RANGE TIME+MAX_H ***;
33000	
33100			next:- nextev;
33200			IF next =/= NONE THEN
33300			WHILE
33400			(IF next == NONE THEN FALSE ELSE
33500			next IN nondiscrete AND next.evtime < time + max_h) DO
33600			next:- next.nextev;
33700	
33800			COMMENT *** SET LARGE_H TO MAX POSSIBLE ***;
33900	
34000			large_h:= IF
34100			(IF  next =/= NONE THEN next.evtime < time + max_h ELSE FALSE)
34200			THEN next.evtime - time
34300			ELSE max_h;
34400	
34500			COMMENT *** IF H SHORTEND ON LAST STEP BY LARGE_H THEN RESTORE VALUE***;
34600	
34700			IF save_h > 0.0
34800			THEN
34900			BEGIN
35000			    h:= save_h;
35100			    save_h:= -1;
35200			END;
35300	
35400			COMMENT *** TEST H SAVE AND REASSIGN IF REQUIRED ***;
35500	
35600			IF h > large_h
35700			THEN
35800			BEGIN
35900			    save_h:= h;
36000			    h:= large_h;
36100			END;
36200	
36300			COMMENT *** ADVANCE INTEGRATION REDUCE STEP SIZE AS REQUIRED  ***;
36400	
36500			take_a_step(h);
36600			WHILE reduce DO
36700			BEGIN
36800			    h:= h / 2;
36900			    FOR i:= 1 STEP 1 UNTIL n DO state[i]:= initial_state[i];
37000			    take_a_step(h);
37100			END;
37200	
37300			COMMENT *** HAS INTEGRATION CAUSED A STATE EVENT ***;
37400	
37500			IF not_state_event
37600			THEN
37700			BEGIN
37800	
37900			    COMMENT *** NO, SO PROCEDE WITH HOLD ,INCREASE H IF POSSIBLE ***;
38000	
38100			    IF trace THEN trace_values;
38200	
38300			    hold(h);
38400			    t:= t + h;
38500	
38600			    IF error_value < 0.07*error_limit AND h <= max_h/2
38700			    THEN h:= 2*h;
38800			END
38900			ELSE
39000			BEGIN
39100	
39200			    COMMENT *** YES, SO START BINARY SEARCH FOR IT ***;
39300	
39400			    ACTIVATE NEW process DELAY h;
39500			    WHILE h >= 2*min_h DO
39600			    BEGIN
39700				WHILE NOT not_state_event AND h >= 2*min_h DO
39800				BEGIN
39900				    FOR i:= 1 STEP 1 UNTIL n DO state[i]:= initial_state[i];
40000				    h:= h/2;
40100				    take_a_step(h);
40200				END *** WHILE STATE EVENT ***;
40300	
40400				WHILE not_state_event AND h >= 2*min_h DO
40500				BEGIN
40600	
40700				    IF trace THEN trace_values;
40800	
40900				    hold(h);
41000				    FOR i:= 1 STEP 1 UNTIL n DO initial_state[i]:= state[i];
41100				    t:= t + h;   h:= h/2;
41200				    take_a_step(h);
41300				END *** WHILE NOT STATE EVENT ***;
41400	
41500			    END *** WHILE H>=2*MIN H ***;
41600			    IF trace THEN trace_values;
41700			    hold(h);
41800			    t:=t+h;
41900	
42000			    stop:= TRUE;
42100			END *** BINARY SEARCH FOR STATE EVENT ***;
42200	
42300		    END -- WHILE LOOP --;
42400	
42500		END --- INTEGRATE ---;
42600	
42700		REAL PROCEDURE state_(i,act_time);   INTEGER i; REAL act_time;
42800		IF (IF idle OR stop THEN TRUE ELSE evtime - start_time <= t) THEN
42900		state_:= state[i] ELSE
43000		BEGIN   REAL a,b,d,k,reltime;
43100		    reltime:= act_time - start_time;
43200		    d:= evtime - start_time - t;
43300		    k:= (state[i]-initial_state[i])/d;
43400		    a:= ((slope_1[i]-k) + (slope_2[i]-k))/(d*d);
43500		    b:= (slope_2[i]-k)/d - a*(evtime - start_time);
43600		    state_:= (reltime-t)*( (act_time-evtime)*
43700		    (reltime*a+b) + k) + initial_state[i];
43800		END OF STATE_ FOR INTERPOLATION OF STATE VARIABLES;
43900	
44000		REAL PROCEDURE rate_(i,act_time);   INTEGER i; REAL act_time;
44100		IF (IF idle OR stop THEN TRUE ELSE evtime-start_time <= t) THEN
44200		rate_:= rate[i] ELSE
44300		BEGIN   REAL a,b,k;
44400		    k:=(act_time-start_time-t)/(evtime-start_time-t);
44500		    a:=2*(slope_1[i]+slope_2[i]-2*slope_3[i]);
44600		    b:=4*slope_3[i]-3*slope_1[i]-slope_2[i];
44700		    rate_:=(a*k+b)*k+slope_1[i];
44800		END OF RATE_ FOR INTERPOLATION OF RATE VALUES;
44900	
45000		OPTIONS(/P:" --- TAKE A STEP(CONTINUOUS) ---");
45100	
45200		COMMENT----------------------------------------------------------------+
45300		!                                                                      !
45400		!   TAKE A STEP IS THE PROCEDURE WHICH ADVANCES THE INTEGRATION BY     !
45500		!   H. IT ALSO PROVIDES AN ESTIMATE OF THE ERROR IN EACH STATE         !
45600		!   VARIABLE. THE FORM PROVIDED IS THE FOURTH ORDER RUNGE-KUTTA MERSON !
45700		!   METHOD.                                                            !
45800		!                                                                      !
45900		!                             PARAMETERS                               !
46000		!                                                                      !
46100		!   H             :THIS, THE ONLY PARAMETER IS THE STEP SIZE THAT      !
46200		!                  THE INTEGRATION IS TO BE ADVANCED BY.               !
46300		!                                                                      !
46400		!                            VARIABLES                                 !
46500		!                                                                      !
46600		!   H_DIV_2, H_DIV_3, H_DIV_6, H_DIV_8  : THE STEP LENGTH, H, DIVIDED  !
46700		!                                         BY 2, 3, 6, 8 RESPECTIVELY.  !
46800		!   I                                   : FOR LOOP INDEX.              !
46900	        !   SLOPE[1],SLOPE[2] AND SLOPE[3] MUST BE CALCULATED HERE, THEY ARE   !
47000	        !   USED FOR THE INTERPOLATION-PROCEDURES (STATE_[ , ] AND RATE_[ , ]) !
47100		+----------------------------------------------------------------------+
47200	;
47300		PROCEDURE take_a_step(h);
47400		VALUE h;
47500		REAL h;
47600		BEGIN
47700		    REAL h_div_2, h_div_3, h_div_6, h_div_8;
47800		    INTEGER i;
47900	
48000		    h_div_2:= h/2;
48100		    h_div_3:= h/3;
48200		    h_div_6:= h/6;
48300		    h_div_8:= h/8;
48400	
48500		    derivatives(t);
48600		    FOR i:= 1 STEP 1 UNTIL n DO
48700		    BEGIN
48800			slope_1[i]:= rate[i];
48900			state[i]:= initial_state[i] + h_div_3 * slope_1[i];
49000		    END;
49100	
49200		    derivatives(t+h_div_3);
49300		    FOR i:= 1 STEP 1 UNTIL n DO
49400		    BEGIN
49500			slope_2[i]:= rate[i];
49600			state[i]:= initial_state[i] + h_div_6 * (slope_1[i] + slope_2[i]);
49700		    END;
49800	
49900		    derivatives(t+h_div_3);
50000		    FOR i:= 1 STEP 1 UNTIL n DO
50100		    BEGIN
50200			slope_2[i]:= rate[i];
50300			state[i]:= initial_state[i] +  h_div_8 * (slope_1[i] + 3 * slope_2[i]);
50400		    END;
50500	
50600		    derivatives(t+h_div_2);
50700		    FOR i:= 1 STEP 1 UNTIL n DO
50800		    BEGIN
50900			slope_3[i]:= rate[i];
51000			state[i]:= error_vector[i]:= initial_state[i] +
51100			h_div_2 * (slope_1[i] - 3*slope_2[i]+ 4*slope_3[i]);
51200		    END;
51300	
51400		    derivatives(t+h);
51500		    FOR i:= 1 STEP 1 UNTIL n DO
51600		    BEGIN
51700			slope_2[i]:= rate[i];
51800			state[i]:= initial_state[i] +
51900			h_div_6*(slope_1[i] + slope_2[i] + 4 * slope_3[i]);
52000			error_vector[i]:= IF relative
52100			THEN Abs(0.2*(1-state[i]/error_vector[i]))
52200			ELSE Abs(0.2*(error_vector[i]-state[i]));
52300	
52400		    END;
52500	
52600		END --- TAKE A STEP ---;
52700	
52800		OPTIONS(/P:" --- DERIVATIVES(CONTINUOUS) ---");
52900	
53000		COMMENT----------------------------------------------------------------+
53100		!                                                                      !
53200		!   DERIVATIVES IS THE PROCEDURE IN WHICH THE USER SPECIFIES THE       !
53300		!   STATE EQUATIONS WHICH GOVERN THE DIFFERENTIAL DYNAMICS OF THE      !
53400		!   CONTINUOUS OBJECT. THE FORM PROVIDED OUTPUTS AN ERROR MESSAGE      !
53500		!   AND ABORTS THE INTEGRATION.                                        !
53600		!                                                                      !
53700		!                             PARAMETERS                               !
53800		!   T             :THIS IS THE TIME RELATIVE TO THE INTEGRATION.       !
53900		!                                                                      !
54000		+----------------------------------------------------------------------+
54100	;
54200		PROCEDURE derivatives(t);
54300		VALUE t;   REAL t;
54400		BEGIN
54500		    abort("? CADSIM ERROR 003 - No PROC DERIVATIVES"
54600		    " HAS BEEN SPECIFIED");
54700		END --- DERIVATIVES ---;
54800	
54900	
55000		max_h:= 100;   min_h:= 0.01;
55100		INNER;
55200	
55300	    END ---*** CONTINUOUS ***---;
55400	
55500	
55600	END **** CADSIM ****;