|   | 3 | == Console1/Console1/Console1.f90 == | 
                  
                          |   | 4 | {{{ | 
                  
                          |   | 5 | #!fortran | 
                  
                          |   | 6 | real x,dx,s | 
                  
                          |   | 7 | !open(1,file='y1.dat', status='unknown') | 
                  
                          |   | 8 | a=-0.5 | 
                  
                          |   | 9 | b=0.5 | 
                  
                          |   | 10 | n=100 | 
                  
                          |   | 11 | dx=(b-a)/n  | 
                  
                          |   | 12 | S1=0 | 
                  
                          |   | 13 | S2=0 | 
                  
                          |   | 14 | S3=0 | 
                  
                          |   | 15 | do i=0,n | 
                  
                          |   | 16 |        S1=S1+(f(i*dx+a)+f((i+1)*dx+a))/2*dx | 
                  
                          |   | 17 |        S2=S2+ f(i*dx+a)*dx  | 
                  
                          |   | 18 | if(mod(i,2)==0) then  | 
                  
                          |   | 19 |        S3=S3+dx/3*2*f(i*dx+a) | 
                  
                          |   | 20 |    else  | 
                  
                          |   | 21 |        S3=S3+dx/3*4*f(i*dx+a) | 
                  
                          |   | 22 |    end if        | 
                  
                          |   | 23 |               | 
                  
                          |   | 24 | end do  | 
                  
                          |   | 25 | print *, 'Trap = ', S1 , 'Rec=', S2 , 'Sim=', S3 | 
                  
                          |   | 26 | pause | 
                  
                          |   | 27 |  | 
                  
                          |   | 28 | contains | 
                  
                          |   | 29 |  | 
                  
                          |   | 30 | function f(x)  | 
                  
                          |   | 31 | f=1/sqrt(1-x**2) | 
                  
                          |   | 32 | end function f | 
                  
                          |   | 33 | end | 
                  
                          |   | 34 | }}} | 
                  
                          |   | 35 |  | 
                  
                          |   | 36 | == Console10/Console10/Console10.f90 == | 
                  
                          |   | 37 | {{{ | 
                  
                          |   | 38 | #!fortran | 
                  
                          |   | 39 | !  Osc_rk_0.f90  | 
                  
                          |   | 40 | ! | 
                  
                          |   | 41 | !  Lab #3 | 
                  
                          |   | 42 | !  metod Runge-Kuta. Oscillator. | 
                  
                          |   | 43 |  | 
                  
                          |   | 44 | external fct, out | 
                  
                          |   | 45 | real aux(8,4) | 
                  
                          |   | 46 | common /a/  dt, k, n, vm, A, om, rl, x, y, r, u(1), u(2), u(3), u(4) | 
                  
                          |   | 47 |  | 
                  
                          |   | 48 | real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/) | 
                  
                          |   | 49 |  | 
                  
                          |   | 50 |  | 
                  
                          |   | 51 | integer :: nd = 2 | 
                  
                          |   | 52 | real :: dt = 2*3.14159265  | 
                  
                          |   | 53 | real :: t | 
                  
                          |   | 54 | real :: c = 3.0e10 | 
                  
                          |   | 55 | !B0=me*c*om/ez  | 
                  
                          |   | 56 | B0 = 200.0 | 
                  
                          |   | 57 | E0 = 5000 | 
                  
                          |   | 58 | eq = 4.8e-10 | 
                  
                          |   | 59 | me = 9.1e-28 | 
                  
                          |   | 60 | om = eq*B0/me*c | 
                  
                          |   | 61 | A = 2.0 | 
                  
                          |   | 62 | vm = A*om | 
                  
                          |   | 63 | n = 1  | 
                  
                          |   | 64 | k = 1 | 
                  
                          |   | 65 | b = B0/B0 | 
                  
                          |   | 66 | rl = c/om | 
                  
                          |   | 67 | u(1) = vm/vm | 
                  
                          |   | 68 | u(2) = vm/vm | 
                  
                          |   | 69 | u(3) = x/rl | 
                  
                          |   | 70 | u(4) = y/rl | 
                  
                          |   | 71 |  | 
                  
                          |   | 72 |  | 
                  
                          |   | 73 |  | 
                  
                          |   | 74 |   | 
                  
                          |   | 75 | open (unit=2, file = '1_rk.dat', status='unknown')  | 
                  
                          |   | 76 |  | 
                  
                          |   | 77 | write(*,*) 'x=', u(3), 'vx=', u(1), 'y=', u(4), 'vy=', u(2) | 
                  
                          |   | 78 | call rkgs(pr, u, du, nd, ih, fct, out, aux) | 
                  
                          |   | 79 |  | 
                  
                          |   | 80 | write(2,*) 'ih=', ih | 
                  
                          |   | 81 |  | 
                  
                          |   | 82 |  | 
                  
                          |   | 83 | end | 
                  
                          |   | 84 |  | 
                  
                          |   | 85 | subroutine fct (t, u, du) | 
                  
                          |   | 86 |  | 
                  
                          |   | 87 | real u(3), du(3), u(4), du(4) | 
                  
                          |   | 88 | common /a/ dt, k, n, vm, A, om, rl, x, y, r,  u(1), u(2), u(3), u(4) | 
                  
                          |   | 89 |    | 
                  
                          |   | 90 | du(1) = -u(3) | 
                  
                          |   | 91 | du(2) = u(4) | 
                  
                          |   | 92 | du(3) = u(1) | 
                  
                          |   | 93 | du(4) = u(2)               | 
                  
                          |   | 94 |                                  | 
                  
                          |   | 95 | return  | 
                  
                          |   | 96 | end   | 
                  
                          |   | 97 |  | 
                  
                          |   | 98 | subroutine  out (t, u, du, ih, nd, pr) | 
                  
                          |   | 99 |  | 
                  
                          |   | 100 | common /a/ dt, k, n, om, vm, A, rl, x, y, r,  u(1), u(2), u(3), u(4) | 
                  
                          |   | 101 |  | 
                  
                          |   | 102 | if (t.ge.k*(dt/20.)) then | 
                  
                          |   | 103 |  | 
                  
                          |   | 104 |  write (4,1) t, u(1), u(2), u(3), u(4)    | 
                  
                          |   | 105 |   k = k + 1 | 
                  
                          |   | 106 |    else | 
                  
                          |   | 107 |     end if | 
                  
                          |   | 108 |  | 
                  
                          |   | 109 | 1 format (f8.3, f8.3, f8.3) | 
                  
                          |   | 110 |  | 
                  
                          |   | 111 | return | 
                  
                          |   | 112 | end  | 
                  
                          |   | 113 |  | 
                  
                          |   | 114 |  | 
                  
                          |   | 115 |       SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)                | 
                  
                          |   | 116 | !                                                                        | 
                  
                          |   | 117 | !                                                                        | 
                  
                          |   | 118 |       DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)             | 
                  
                          |   | 119 |       DO 1 I=1,NDIM                                                      | 
                  
                          |   | 120 |     1 AUX(8,I)=.06666667*DERY(I)                                         | 
                  
                          |   | 121 |       X=PRMT(1)                                                          | 
                  
                          |   | 122 |       XEND=PRMT(2)                                                       | 
                  
                          |   | 123 |       H=PRMT(3)                                                          | 
                  
                          |   | 124 |       PRMT(5)=0.                                                         | 
                  
                          |   | 125 |       CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 126 | !                                                                        | 
                  
                          |   | 127 | !     ERROR TEST                                                         | 
                  
                          |   | 128 |       IF(H*(XEND-X))38,37,2                                              | 
                  
                          |   | 129 | !                                                                        | 
                  
                          |   | 130 | !     PREPARATIONS FOR RUNGE-KUTTA METHOD                                | 
                  
                          |   | 131 |     2 A(1)=.5                                                            | 
                  
                          |   | 132 |       A(2)=.2928932                                                      | 
                  
                          |   | 133 |       A(3)=1.707107                                                      | 
                  
                          |   | 134 |       A(4)=.1666667                                                      | 
                  
                          |   | 135 |       B(1)=2.                                                            | 
                  
                          |   | 136 |       B(2)=1.                                                            | 
                  
                          |   | 137 |       B(3)=1.                                                            | 
                  
                          |   | 138 |       B(4)=2.                                                            | 
                  
                          |   | 139 |       C(1)=.5                                                            | 
                  
                          |   | 140 |       C(2)=.2928932                                                      | 
                  
                          |   | 141 |       C(3)=1.707107                                                      | 
                  
                          |   | 142 |       C(4)=.5                                                            | 
                  
                          |   | 143 | !                                                                        | 
                  
                          |   | 144 | !     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                             | 
                  
                          |   | 145 |       DO 3 I=1,NDIM                                                      | 
                  
                          |   | 146 |       AUX(1,I)=Y(I)                                                      | 
                  
                          |   | 147 |       AUX(2,I)=DERY(I)                                                   | 
                  
                          |   | 148 |       AUX(3,I)=0.                                                        | 
                  
                          |   | 149 |     3 AUX(6,I)=0.                                                        | 
                  
                          |   | 150 |       IREC=0                                                             | 
                  
                          |   | 151 |       H=H+H                                                              | 
                  
                          |   | 152 |       IHLF=-1                                                            | 
                  
                          |   | 153 |       ISTEP=0                                                            | 
                  
                          |   | 154 |       IEND=0                                                             | 
                  
                          |   | 155 | !                                                                       | 
                  
                          |   | 156 | !                                                                    | 
                  
                          |   | 157 | !     START OF A RUNGE-KUTTA STEP                                        | 
                  
                          |   | 158 |     4 IF((X+H-XEND)*H)7,6,5                                              | 
                  
                          |   | 159 |     5 H=XEND-X                                                           | 
                  
                          |   | 160 |     6 IEND=1                                                             | 
                  
                          |   | 161 | !                                                                       | 
                  
                          |   | 162 | !     RECORDING OF INITIAL VALUES OF THIS STEP                           | 
                  
                          |   | 163 |     7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                 | 
                  
                          |   | 164 |       IF(PRMT(5))40,8,40                                                 | 
                  
                          |   | 165 |     8 ITEST=0                                                            | 
                  
                          |   | 166 |     9 ISTEP=ISTEP+1                                                      | 
                  
                          |   | 167 | !     START OF INNERMOST RUNGE-KUTTA LOOP                                | 
                  
                          |   | 168 |       J=1                                                                | 
                  
                          |   | 169 |    10 AJ=A(J)                                                            | 
                  
                          |   | 170 |       BJ=B(J)                                                            | 
                  
                          |   | 171 |       CJ=C(J)                                                            | 
                  
                          |   | 172 |       DO 11 I=1,NDIM                                                     | 
                  
                          |   | 173 |       R1=H*DERY(I)                                                       | 
                  
                          |   | 174 |       R2=AJ*(R1-BJ*AUX(6,I))                                             | 
                  
                          |   | 175 |       Y(I)=Y(I)+R2                                                       | 
                  
                          |   | 176 |       R2=R2+R2+R2                                                        | 
                  
                          |   | 177 |    11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                         | 
                  
                          |   | 178 |       IF(J-4)12,15,15                                                    | 
                  
                          |   | 179 |    12 J=J+1                                                              | 
                  
                          |   | 180 |       IF(J-3)13,14,13                                                    | 
                  
                          |   | 181 |    13 X=X+.5*H                                                           | 
                  
                          |   | 182 |    14 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 183 |       GOTO 10                                                            | 
                  
                          |   | 184 | !     END OF INNERMOST RUNGE-KUTTA LOOP                                  | 
                  
                          |   | 185 | !     TEST OF ACCURACY                                                   | 
                  
                          |   | 186 |    15 IF(ITEST)16,16,20                                                  | 
                  
                          |   | 187 | !                                                                        | 
                  
                          |   | 188 | !     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY    | 
                  
                          |   | 189 |    16 DO 17 I=1,NDIM                                                     | 
                  
                          |   | 190 |    17 AUX(4,I)=Y(I)                                                      | 
                  
                          |   | 191 |       ITEST=1                                                            | 
                  
                          |   | 192 |       ISTEP=ISTEP+ISTEP-2                                                | 
                  
                          |   | 193 |    18 IHLF=IHLF+1 | 
                  
                          |   | 194 |       X=X-H                                                              | 
                  
                          |   | 195 |       H=.5*H                                                             | 
                  
                          |   | 196 |       DO 19 I=1,NDIM                                                     | 
                  
                          |   | 197 |       Y(I)=AUX(1,I)                                                      | 
                  
                          |   | 198 |       DERY(I)=AUX(2,I)                                                   | 
                  
                          |   | 199 |    19 AUX(6,I)=AUX(3,I)                                                  | 
                  
                          |   | 200 |       GOTO 9                                                             | 
                  
                          |   | 201 | !                                                                        | 
                  
                          |   | 202 | !     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                    | 
                  
                          |   | 203 |    20 IMOD=ISTEP/2                                                       | 
                  
                          |   | 204 |       IF(ISTEP-IMOD-IMOD)21,23,21                                        | 
                  
                          |   | 205 |    21 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 206 |       DO 22 I=1,NDIM                                                     | 
                  
                          |   | 207 |       AUX(5,I)=Y(I)                                                      | 
                  
                          |   | 208 |    22 AUX(7,I)=DERY(I)                                                   | 
                  
                          |   | 209 |       GOTO 9                                                             | 
                  
                          |   | 210 | !                                                                        | 
                  
                          |   | 211 | !     COMPUTATION OF TEST VALUE DELT                                     | 
                  
                          |   | 212 |    23 DELT=0.                                                            | 
                  
                          |   | 213 |       DO 24 I=1,NDIM                                                     | 
                  
                          |   | 214 |    24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                              | 
                  
                          |   | 215 |       IF(DELT-PRMT(4))28,28,25                                           | 
                  
                          |   | 216 | !                                                                        | 
                  
                          |   | 217 | !     ERROR IS TOO GREAT                                                 | 
                  
                          |   | 218 |    25 IF(IHLF-10)26,36,36                                                | 
                  
                          |   | 219 |    26 DO 27 I=1,NDIM                                                     | 
                  
                          |   | 220 |    27 AUX(4,I)=AUX(5,I)                                                  | 
                  
                          |   | 221 |       ISTEP=ISTEP+ISTEP-4                                                | 
                  
                          |   | 222 |       X=X-H                                                              | 
                  
                          |   | 223 |       IEND=0                                                             | 
                  
                          |   | 224 |       GOTO 18                                                            | 
                  
                          |   | 225 | !                                                                        | 
                  
                          |   | 226 | !     RESULT VALUES ARE GOOD                                             | 
                  
                          |   | 227 |    28 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 228 |       DO 29 I=1,NDIM                                                     | 
                  
                          |   | 229 |       AUX(1,I)=Y(I)                                                      | 
                  
                          |   | 230 |       AUX(2,I)=DERY(I)                                                   | 
                  
                          |   | 231 |       AUX(3,I)=AUX(6,I)                                                  | 
                  
                          |   | 232 |       Y(I)=AUX(5,I)                                                      | 
                  
                          |   | 233 |    29 DERY(I)=AUX(7,I)                                                   | 
                  
                          |   | 234 |       CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                               | 
                  
                          |   | 235 |       IF(PRMT(5))40,30,40                                                | 
                  
                          |   | 236 |    30 DO 31 I=1,NDIM                                                     | 
                  
                          |   | 237 |       Y(I)=AUX(1,I)                                                      | 
                  
                          |   | 238 |    31 DERY(I)=AUX(2,I)                                                   | 
                  
                          |   | 239 |       IREC=IHLF                                                          | 
                  
                          |   | 240 |       IF(IEND)32,32,39                                                   | 
                  
                          |   | 241 | !                                                                        | 
                  
                          |   | 242 | !     INCREMENT GETS DOUBLED                                             | 
                  
                          |   | 243 |    32 IHLF=IHLF-1                                                        | 
                  
                          |   | 244 |       ISTEP=ISTEP/2                                                      | 
                  
                          |   | 245 |       H=H+H                                                              | 
                  
                          |   | 246 |       IF(IHLF)4,33,33                                                    | 
                  
                          |   | 247 |    33 IMOD=ISTEP/2                                                       | 
                  
                          |   | 248 |       IF(ISTEP-IMOD-IMOD)4,34,4                                          | 
                  
                          |   | 249 |    34 IF(DELT-.02*PRMT(4))35,35,4                                        | 
                  
                          |   | 250 |    35 IHLF=IHLF-1                                                        | 
                  
                          |   | 251 |       ISTEP=ISTEP/2                                                      | 
                  
                          |   | 252 |       H=H+H                                                              | 
                  
                          |   | 253 |       GOTO 4                                                             | 
                  
                          |   | 254 | !                                                                        | 
                  
                          |   | 255 | !     RETURNS TO CALLING PROGRAM                                         | 
                  
                          |   | 256 |    36 IHLF=11                                                            | 
                  
                          |   | 257 |       CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 258 |       GOTO 39                                                            | 
                  
                          |   | 259 |    37 IHLF=12                                                            | 
                  
                          |   | 260 |       GOTO 39                                                            | 
                  
                          |   | 261 |    38 IHLF=13                                                            | 
                  
                          |   | 262 |    39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                 | 
                  
                          |   | 263 |    40 RETURN                                                             | 
                  
                          |   | 264 |       END                                                                | 
                  
                          |   | 265 |  | 
                  
                          |   | 266 |  | 
                  
                          |   | 267 | }}} | 
                  
                          |   | 268 |  | 
                  
                          |   | 269 | == Console11/Console11/Console11.f90 == | 
                  
                          |   | 270 | {{{ | 
                  
                          |   | 271 | #!fortran | 
                  
                          |   | 272 | !  Console11.f90  | 
                  
                          |   | 273 | ! | 
                  
                          |   | 274 | !  FUNCTIONS: | 
                  
                          |   | 275 | !  Console11 - Entry point of console application. | 
                  
                          |   | 276 | ! | 
                  
                          |   | 277 |  | 
                  
                          |   | 278 | !**************************************************************************** | 
                  
                          |   | 279 | ! | 
                  
                          |   | 280 | !  PROGRAM: Console11 | 
                  
                          |   | 281 | ! | 
                  
                          |   | 282 | !  PURPOSE:  Entry point for the console application. | 
                  
                          |   | 283 | ! | 
                  
                          |   | 284 | !**************************************************************************** | 
                  
                          |   | 285 |  | 
                  
                          |   | 286 |     program Console11 | 
                  
                          |   | 287 |  | 
                  
                          |   | 288 |     implicit none | 
                  
                          |   | 289 |  | 
                  
                          |   | 290 |     ! Variables | 
                  
                          |   | 291 |  | 
                  
                          |   | 292 |     ! Body of Console11 | 
                  
                          |   | 293 |     print *, 'Hello World' | 
                  
                          |   | 294 |  | 
                  
                          |   | 295 |     end program Console11 | 
                  
                          |   | 296 |  | 
                  
                          |   | 297 | }}} | 
                  
                          |   | 298 |  | 
                  
                          |   | 299 | == Console12/Console12/Console12.f90 == | 
                  
                          |   | 300 | {{{ | 
                  
                          |   | 301 | #!fortran | 
                  
                          |   | 302 | !  Console12.f90  | 
                  
                          |   | 303 | ! | 
                  
                          |   | 304 | !  FUNCTIONS: | 
                  
                          |   | 305 | !  Console12 - Entry point of console application. | 
                  
                          |   | 306 | ! | 
                  
                          |   | 307 |  | 
                  
                          |   | 308 | !**************************************************************************** | 
                  
                          |   | 309 | ! | 
                  
                          |   | 310 | !  PROGRAM: Console12 | 
                  
                          |   | 311 | ! | 
                  
                          |   | 312 | !  PURPOSE:  Entry point for the console application. | 
                  
                          |   | 313 | ! | 
                  
                          |   | 314 | !**************************************************************************** | 
                  
                          |   | 315 |  | 
                  
                          |   | 316 |     program Console12 | 
                  
                          |   | 317 |  | 
                  
                          |   | 318 |     implicit none | 
                  
                          |   | 319 |  | 
                  
                          |   | 320 |     ! Variables | 
                  
                          |   | 321 |  | 
                  
                          |   | 322 |     ! Body of Console12 | 
                  
                          |   | 323 |     print *, 'Hello World' | 
                  
                          |   | 324 |  | 
                  
                          |   | 325 |     end program Console12 | 
                  
                          |   | 326 |  | 
                  
                          |   | 327 | }}} | 
                  
                          |   | 328 |  | 
                  
                          |   | 329 | == Console13/Console13/Console13.f90 == | 
                  
                          |   | 330 | {{{ | 
                  
                          |   | 331 | #!fortran | 
                  
                          |   | 332 | integer :: n=1000 | 
                  
                          |   | 333 | real(8) :: dt, t, vm, b, rl, om, x, y, vx, vy, B0 | 
                  
                          |   | 334 | real(8) :: eq=-4.8e-10, me=9.1e-28, c = 3.0e10 | 
                  
                          |   | 335 | !B0=me*c*om/ez  | 
                  
                          |   | 336 | B0 = 200.0 | 
                  
                          |   | 337 | !E0 = 5000 | 
                  
                          |   | 338 | !A = 1e-30 | 
                  
                          |   | 339 | x=1.0 | 
                  
                          |   | 340 | y=1.0 | 
                  
                          |   | 341 | vx=1000.0 | 
                  
                          |   | 342 | vy=1000.0 | 
                  
                          |   | 343 | a=100 | 
                  
                          |   | 344 | !b = B0/B0 | 
                  
                          |   | 345 | !om = eq*B0/me*c | 
                  
                          |   | 346 | !vm = A*om | 
                  
                          |   | 347 | !print *, vm | 
                  
                          |   | 348 | !rl = c/om | 
                  
                          |   | 349 | !vx = vm/vm | 
                  
                          |   | 350 | !vy = vm/vm | 
                  
                          |   | 351 | !x = x/rl | 
                  
                          |   | 352 | !y = y/rl | 
                  
                          |   | 353 | !print *, rl | 
                  
                          |   | 354 | t = 1.0e-8 | 
                  
                          |   | 355 | dt = t/float(n) | 
                  
                          |   | 356 |  | 
                  
                          |   | 357 | open (1,file='1',status='unknown') | 
                  
                          |   | 358 |  | 
                  
                          |   | 359 | do i=1,n  | 
                  
                          |   | 360 |      vx = vx+(B0*eq)/(c*me)*vy*dt | 
                  
                          |   | 361 |      vy = vy-(B0*eq)/(c*me)*vx*dt | 
                  
                          |   | 362 |      x = x+vx*dt | 
                  
                          |   | 363 |      y = y+vy*dt | 
                  
                          |   | 364 |    write(1,*) x, vx | 
                  
                          |   | 365 | end do | 
                  
                          |   | 366 |  | 
                  
                          |   | 367 | pause | 
                  
                          |   | 368 | end program | 
                  
                          |   | 369 | }}} | 
                  
                          |   | 370 |  | 
                  
                          |   | 371 | == Console15/Console15/Console15.f90 == | 
                  
                          |   | 372 | {{{ | 
                  
                          |   | 373 | #!fortran | 
                  
                          |   | 374 | x = (/1,5/) | 
                  
                          |   | 375 | print *, x | 
                  
                          |   | 376 | end | 
                  
                          |   | 377 | }}} | 
                  
                          |   | 378 |  | 
                  
                          |   | 379 | == Console16/Console16/Console16.f90 == | 
                  
                          |   | 380 | {{{ | 
                  
                          |   | 381 | #!fortran | 
                  
                          |   | 382 | integer :: n=1000 | 
                  
                          |   | 383 | real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b | 
                  
                          |   | 384 | real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0  | 
                  
                          |   | 385 | x = 1.0 | 
                  
                          |   | 386 | y = 1.0 | 
                  
                          |   | 387 | a = 0.5 | 
                  
                          |   | 388 |  | 
                  
                          |   | 389 | vx = 1.0e09 | 
                  
                          |   | 390 | vy = 1.0e09 | 
                  
                          |   | 391 | E0=0.0 | 
                  
                          |   | 392 | B0 = 200.0 | 
                  
                          |   | 393 |  | 
                  
                          |   | 394 | b=(2*3.14*6)/(1.0e09) | 
                  
                          |   | 395 |  | 
                  
                          |   | 396 | !B0 = 1e-25 | 
                  
                          |   | 397 | !E0 = 5000 | 
                  
                          |   | 398 | !A = 1e-30 | 
                  
                          |   | 399 | !b = B0/B0 | 
                  
                          |   | 400 | !om = eq*B0/me*c | 
                  
                          |   | 401 | !vm = A*om | 
                  
                          |   | 402 | !print *, vm | 
                  
                          |   | 403 | !rl = c/om | 
                  
                          |   | 404 | !vx = vm/vm | 
                  
                          |   | 405 | !vy = vm/vm | 
                  
                          |   | 406 | !x = x/rl | 
                  
                          |   | 407 | !y = y/rl | 
                  
                          |   | 408 | !print *, rl | 
                  
                          |   | 409 | t = 1.00e-08 | 
                  
                          |   | 410 | dt = t/float(n) | 
                  
                          |   | 411 |  | 
                  
                          |   | 412 |  | 
                  
                          |   | 413 | open (1,file='1',status='unknown') | 
                  
                          |   | 414 |  | 
                  
                          |   | 415 | do i=1,n | 
                  
                          |   | 416 |     | 
                  
                          |   | 417 |    Bz = B0*(1+a*y)  | 
                  
                          |   | 418 |    !Bz = B0 | 
                  
                          |   | 419 |    x = x + vx*dt | 
                  
                          |   | 420 |    y = y + vy*dt | 
                  
                          |   | 421 |  | 
                  
                          |   | 422 |    !vx = vx + (vy*eq*Bz)/(c*me)*dt | 
                  
                          |   | 423 |    !vy = vy - (vx*eq*Bz)/(c*me)*dt | 
                  
                          |   | 424 |  | 
                  
                          |   | 425 |     | 
                  
                          |   | 426 |    if(i*dt.gt.b) then  | 
                  
                          |   | 427 | E0 = 5.0 | 
                  
                          |   | 428 | end if | 
                  
                          |   | 429 |     | 
                  
                          |   | 430 |   vx = vx +(eq*E0)/(me)*dt+(vy*eq*Bz)/(c*me)*dt | 
                  
                          |   | 431 |    vy = vy + (eq*E0)/(me)*dt-(vx*eq*Bz)/(c*me)*dt | 
                  
                          |   | 432 |   ! x = x + vx*dt | 
                  
                          |   | 433 |    !y = y + vy*dt | 
                  
                          |   | 434 |  | 
                  
                          |   | 435 |    !print *, vx, vy | 
                  
                          |   | 436 |  | 
                  
                          |   | 437 |  | 
                  
                          |   | 438 |  | 
                  
                          |   | 439 |    | 
                  
                          |   | 440 |  | 
                  
                          |   | 441 |           write (1,*) x, y  | 
                  
                          |   | 442 |  | 
                  
                          |   | 443 |     end do | 
                  
                          |   | 444 |  pause   | 
                  
                          |   | 445 | end program | 
                  
                          |   | 446 | }}} | 
                  
                          |   | 447 |  | 
                  
                          |   | 448 | == Console17/Console17/Console17.f90 == | 
                  
                          |   | 449 | {{{ | 
                  
                          |   | 450 | #!fortran | 
                  
                          |   | 451 | integer :: N = 1500 | 
                  
                          |   | 452 | real(8) :: x = 2.5, y = 0.5, B0 = 500.0, We = 4.0054e-08, b = 2.0, pi = 3.14159, m = 9.1e-28, e = -4.8e-10, c = 2.99e10 | 
                  
                          |   | 453 | real(8) vx, vy, v, dt, t | 
                  
                          |   | 454 |  | 
                  
                          |   | 455 | v = sqrt(2*We/m) | 
                  
                          |   | 456 | vx = sin(pi/4.0)*v | 
                  
                          |   | 457 | vy = cos(pi/4.0)*v | 
                  
                          |   | 458 | t = 5e-9 | 
                  
                          |   | 459 | dt=t/float(N) | 
                  
                          |   | 460 | !Norma | 
                  
                          |   | 461 | !tau = t*10e8 | 
                  
                          |   | 462 | !me = m*1e28 | 
                  
                          |   | 463 | !qe = e*1e10 | 
                  
                          |   | 464 | !vl = c*1e-10 | 
                  
                          |   | 465 | !ve = v/c | 
                  
                          |   | 466 | !vx = sin(pi/4.0)*v | 
                  
                          |   | 467 | !vy = cos(pi/4.0)*v | 
                  
                          |   | 468 | !dt = tau/float(N) | 
                  
                          |   | 469 |  | 
                  
                          |   | 470 | open (1,file='mov',status='replace') | 
                  
                          |   | 471 | do i=1,N | 
                  
                          |   | 472 |    Bz = B0*(1-b**2*(x**2+y**2)) | 
                  
                          |   | 473 |    x = x + vx*dt | 
                  
                          |   | 474 |    y = y + vy*dt | 
                  
                          |   | 475 |    vx = vx + vy/m*Bz*e/c*dt | 
                  
                          |   | 476 |    vy = vy - vx/m*Bz*e/c*dt | 
                  
                          |   | 477 |    print *, x, y | 
                  
                          |   | 478 |    pause | 
                  
                          |   | 479 |    write(1,*) y, x | 
                  
                          |   | 480 | end do | 
                  
                          |   | 481 | end program | 
                  
                          |   | 482 | }}} | 
                  
                          |   | 483 |  | 
                  
                          |   | 484 | == Console18/Console18/Console18.f90 == | 
                  
                          |   | 485 | {{{ | 
                  
                          |   | 486 | #!fortran | 
                  
                          |   | 487 | !  Osc_rk_0.f90  | 
                  
                          |   | 488 | ! | 
                  
                          |   | 489 | !  Lab #3 | 
                  
                          |   | 490 | !  metod Runge-Kuta. Oscillator. | 
                  
                          |   | 491 |  | 
                  
                          |   | 492 | external fct, out | 
                  
                          |   | 493 | real aux(8,2) | 
                  
                          |   | 494 | common /a/ b, lam, f0, dt, k, n | 
                  
                          |   | 495 |  | 
                  
                          |   | 496 | real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/) | 
                  
                          |   | 497 | real :: u(2) = (/0., 0.5/)  | 
                  
                          |   | 498 | real :: du(2) = (/0.5, 0.5/) | 
                  
                          |   | 499 | real :: b=0.003 | 
                  
                          |   | 500 | integer :: nd = 2 | 
                  
                          |   | 501 |            | 
                  
                          |   | 502 | real :: dt = 2*3.14159265  | 
                  
                          |   | 503 |     | 
                  
                          |   | 504 | ! real :: lam = 0.1 | 
                  
                          |   | 505 | real :: f0 = 0.00 | 
                  
                          |   | 506 |  | 
                  
                          |   | 507 | n = 1  | 
                  
                          |   | 508 | k = 1 | 
                  
                          |   | 509 |   | 
                  
                          |   | 510 | open (unit=2, file = 'osc_rk.dat', status='new')  | 
                  
                          |   | 511 |  | 
                  
                          |   | 512 | write(*,*) 'x=', u(1), 'v=', u(2) | 
                  
                          |   | 513 | pause | 
                  
                          |   | 514 | call rkgs(pr, u, du, nd, ih, fct, out, aux) | 
                  
                          |   | 515 |  | 
                  
                          |   | 516 | write(2,*) 'ih=', ih | 
                  
                          |   | 517 |  | 
                  
                          |   | 518 | stop | 
                  
                          |   | 519 | end | 
                  
                          |   | 520 |  | 
                  
                          |   | 521 | subroutine fct (t, u, du) | 
                  
                          |   | 522 |  | 
                  
                          |   | 523 | real u(2), du(2), lam  | 
                  
                          |   | 524 | common /a/ b, lam, f0, dt, k, n | 
                  
                          |   | 525 |  | 
                  
                          |   | 526 | du(1) = -u(2)-2.*b*u(1)                   | 
                  
                          |   | 527 | du(2) = u(1) | 
                  
                          |   | 528 |                                  | 
                  
                          |   | 529 | return  | 
                  
                          |   | 530 | end   | 
                  
                          |   | 531 |  | 
                  
                          |   | 532 | subroutine  out (t, u, du, ih, nd, pr) | 
                  
                          |   | 533 | real u(2), du(2), pr(5), lam | 
                  
                          |   | 534 | common /a/ b, lam, f0, dt, k, n | 
                  
                          |   | 535 |  | 
                  
                          |   | 536 | if (t.ge.k*(dt/20.)) then | 
                  
                          |   | 537 |  | 
                  
                          |   | 538 |  write (2,1) t,u(1),u(2)  !  | 
                  
                          |   | 539 |   k = k + 1 | 
                  
                          |   | 540 |    else | 
                  
                          |   | 541 |     end if | 
                  
                          |   | 542 |  | 
                  
                          |   | 543 | 1 format (f8.3, f8.3, f8.3) | 
                  
                          |   | 544 |  | 
                  
                          |   | 545 | return | 
                  
                          |   | 546 | end  | 
                  
                          |   | 547 |  | 
                  
                          |   | 548 |  | 
                  
                          |   | 549 |       SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)                | 
                  
                          |   | 550 | !                                                                        | 
                  
                          |   | 551 | !                                                                        | 
                  
                          |   | 552 |       DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)             | 
                  
                          |   | 553 |       DO 1 I=1,NDIM                                                      | 
                  
                          |   | 554 |     1 AUX(8,I)=.06666667*DERY(I)                                         | 
                  
                          |   | 555 |       X=PRMT(1)                                                          | 
                  
                          |   | 556 |       XEND=PRMT(2)                                                       | 
                  
                          |   | 557 |       H=PRMT(3)                                                          | 
                  
                          |   | 558 |       PRMT(5)=0.                                                         | 
                  
                          |   | 559 |       CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 560 | !                                                                        | 
                  
                          |   | 561 | !     ERROR TEST                                                         | 
                  
                          |   | 562 |       IF(H*(XEND-X))38,37,2                                              | 
                  
                          |   | 563 | !                                                                        | 
                  
                          |   | 564 | !     PREPARATIONS FOR RUNGE-KUTTA METHOD                                | 
                  
                          |   | 565 |     2 A(1)=.5                                                            | 
                  
                          |   | 566 |       A(2)=.2928932                                                      | 
                  
                          |   | 567 |       A(3)=1.707107                                                      | 
                  
                          |   | 568 |       A(4)=.1666667                                                      | 
                  
                          |   | 569 |       B(1)=2.                                                            | 
                  
                          |   | 570 |       B(2)=1.                                                            | 
                  
                          |   | 571 |       B(3)=1.                                                            | 
                  
                          |   | 572 |       B(4)=2.                                                            | 
                  
                          |   | 573 |       C(1)=.5                                                            | 
                  
                          |   | 574 |       C(2)=.2928932                                                      | 
                  
                          |   | 575 |       C(3)=1.707107                                                      | 
                  
                          |   | 576 |       C(4)=.5                                                            | 
                  
                          |   | 577 | !                                                                        | 
                  
                          |   | 578 | !     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                             | 
                  
                          |   | 579 |       DO 3 I=1,NDIM                                                      | 
                  
                          |   | 580 |       AUX(1,I)=Y(I)                                                      | 
                  
                          |   | 581 |       AUX(2,I)=DERY(I)                                                   | 
                  
                          |   | 582 |       AUX(3,I)=0.                                                        | 
                  
                          |   | 583 |     3 AUX(6,I)=0.                                                        | 
                  
                          |   | 584 |       IREC=0                                                             | 
                  
                          |   | 585 |       H=H+H                                                              | 
                  
                          |   | 586 |       IHLF=-1                                                            | 
                  
                          |   | 587 |       ISTEP=0                                                            | 
                  
                          |   | 588 |       IEND=0                                                             | 
                  
                          |   | 589 | !                                                                       | 
                  
                          |   | 590 | !                                                                    | 
                  
                          |   | 591 | !     START OF A RUNGE-KUTTA STEP                                        | 
                  
                          |   | 592 |     4 IF((X+H-XEND)*H)7,6,5                                              | 
                  
                          |   | 593 |     5 H=XEND-X                                                           | 
                  
                          |   | 594 |     6 IEND=1                                                             | 
                  
                          |   | 595 | !                                                                       | 
                  
                          |   | 596 | !     RECORDING OF INITIAL VALUES OF THIS STEP                           | 
                  
                          |   | 597 |     7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                 | 
                  
                          |   | 598 |       IF(PRMT(5))40,8,40                                                 | 
                  
                          |   | 599 |     8 ITEST=0                                                            | 
                  
                          |   | 600 |     9 ISTEP=ISTEP+1                                                      | 
                  
                          |   | 601 | !     START OF INNERMOST RUNGE-KUTTA LOOP                                | 
                  
                          |   | 602 |       J=1                                                                | 
                  
                          |   | 603 |    10 AJ=A(J)                                                            | 
                  
                          |   | 604 |       BJ=B(J)                                                            | 
                  
                          |   | 605 |       CJ=C(J)                                                            | 
                  
                          |   | 606 |       DO 11 I=1,NDIM                                                     | 
                  
                          |   | 607 |       R1=H*DERY(I)                                                       | 
                  
                          |   | 608 |       R2=AJ*(R1-BJ*AUX(6,I))                                             | 
                  
                          |   | 609 |       Y(I)=Y(I)+R2                                                       | 
                  
                          |   | 610 |       R2=R2+R2+R2                                                        | 
                  
                          |   | 611 |    11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                         | 
                  
                          |   | 612 |       IF(J-4)12,15,15                                                    | 
                  
                          |   | 613 |    12 J=J+1                                                              | 
                  
                          |   | 614 |       IF(J-3)13,14,13                                                    | 
                  
                          |   | 615 |    13 X=X+.5*H                                                           | 
                  
                          |   | 616 |    14 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 617 |       GOTO 10                                                            | 
                  
                          |   | 618 | !     END OF INNERMOST RUNGE-KUTTA LOOP                                  | 
                  
                          |   | 619 | !     TEST OF ACCURACY                                                   | 
                  
                          |   | 620 |    15 IF(ITEST)16,16,20                                                  | 
                  
                          |   | 621 | !                                                                        | 
                  
                          |   | 622 | !     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY    | 
                  
                          |   | 623 |    16 DO 17 I=1,NDIM                                                     | 
                  
                          |   | 624 |    17 AUX(4,I)=Y(I)                                                      | 
                  
                          |   | 625 |       ITEST=1                                                            | 
                  
                          |   | 626 |       ISTEP=ISTEP+ISTEP-2                                                | 
                  
                          |   | 627 |    18 IHLF=IHLF+1 | 
                  
                          |   | 628 |       X=X-H                                                              | 
                  
                          |   | 629 |       H=.5*H                                                             | 
                  
                          |   | 630 |       DO 19 I=1,NDIM                                                     | 
                  
                          |   | 631 |       Y(I)=AUX(1,I)                                                      | 
                  
                          |   | 632 |       DERY(I)=AUX(2,I)                                                   | 
                  
                          |   | 633 |    19 AUX(6,I)=AUX(3,I)                                                  | 
                  
                          |   | 634 |       GOTO 9                                                             | 
                  
                          |   | 635 | !                                                                        | 
                  
                          |   | 636 | !     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                    | 
                  
                          |   | 637 |    20 IMOD=ISTEP/2                                                       | 
                  
                          |   | 638 |       IF(ISTEP-IMOD-IMOD)21,23,21                                        | 
                  
                          |   | 639 |    21 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 640 |       DO 22 I=1,NDIM                                                     | 
                  
                          |   | 641 |       AUX(5,I)=Y(I)                                                      | 
                  
                          |   | 642 |    22 AUX(7,I)=DERY(I)                                                   | 
                  
                          |   | 643 |       GOTO 9                                                             | 
                  
                          |   | 644 | !                                                                        | 
                  
                          |   | 645 | !     COMPUTATION OF TEST VALUE DELT                                     | 
                  
                          |   | 646 |    23 DELT=0.                                                            | 
                  
                          |   | 647 |       DO 24 I=1,NDIM                                                     | 
                  
                          |   | 648 |    24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                              | 
                  
                          |   | 649 |       IF(DELT-PRMT(4))28,28,25                                           | 
                  
                          |   | 650 | !                                                                        | 
                  
                          |   | 651 | !     ERROR IS TOO GREAT                                                 | 
                  
                          |   | 652 |    25 IF(IHLF-10)26,36,36                                                | 
                  
                          |   | 653 |    26 DO 27 I=1,NDIM                                                     | 
                  
                          |   | 654 |    27 AUX(4,I)=AUX(5,I)                                                  | 
                  
                          |   | 655 |       ISTEP=ISTEP+ISTEP-4                                                | 
                  
                          |   | 656 |       X=X-H                                                              | 
                  
                          |   | 657 |       IEND=0                                                             | 
                  
                          |   | 658 |       GOTO 18                                                            | 
                  
                          |   | 659 | !                                                                        | 
                  
                          |   | 660 | !     RESULT VALUES ARE GOOD                                             | 
                  
                          |   | 661 |    28 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 662 |       DO 29 I=1,NDIM                                                     | 
                  
                          |   | 663 |       AUX(1,I)=Y(I)                                                      | 
                  
                          |   | 664 |       AUX(2,I)=DERY(I)                                                   | 
                  
                          |   | 665 |       AUX(3,I)=AUX(6,I)                                                  | 
                  
                          |   | 666 |       Y(I)=AUX(5,I)                                                      | 
                  
                          |   | 667 |    29 DERY(I)=AUX(7,I)                                                   | 
                  
                          |   | 668 |       CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                               | 
                  
                          |   | 669 |       IF(PRMT(5))40,30,40                                                | 
                  
                          |   | 670 |    30 DO 31 I=1,NDIM                                                     | 
                  
                          |   | 671 |       Y(I)=AUX(1,I)                                                      | 
                  
                          |   | 672 |    31 DERY(I)=AUX(2,I)                                                   | 
                  
                          |   | 673 |       IREC=IHLF                                                          | 
                  
                          |   | 674 |       IF(IEND)32,32,39                                                   | 
                  
                          |   | 675 | !                                                                        | 
                  
                          |   | 676 | !     INCREMENT GETS DOUBLED                                             | 
                  
                          |   | 677 |    32 IHLF=IHLF-1                                                        | 
                  
                          |   | 678 |       ISTEP=ISTEP/2                                                      | 
                  
                          |   | 679 |       H=H+H                                                              | 
                  
                          |   | 680 |       IF(IHLF)4,33,33                                                    | 
                  
                          |   | 681 |    33 IMOD=ISTEP/2                                                       | 
                  
                          |   | 682 |       IF(ISTEP-IMOD-IMOD)4,34,4                                          | 
                  
                          |   | 683 |    34 IF(DELT-.02*PRMT(4))35,35,4                                        | 
                  
                          |   | 684 |    35 IHLF=IHLF-1                                                        | 
                  
                          |   | 685 |       ISTEP=ISTEP/2                                                      | 
                  
                          |   | 686 |       H=H+H                                                              | 
                  
                          |   | 687 |       GOTO 4                                                             | 
                  
                          |   | 688 | !                                                                        | 
                  
                          |   | 689 | !     RETURNS TO CALLING PROGRAM                                         | 
                  
                          |   | 690 |    36 IHLF=11                                                            | 
                  
                          |   | 691 |       CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 692 |       GOTO 39                                                            | 
                  
                          |   | 693 |    37 IHLF=12                                                            | 
                  
                          |   | 694 |       GOTO 39                                                            | 
                  
                          |   | 695 |    38 IHLF=13                                                            | 
                  
                          |   | 696 |    39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                 | 
                  
                          |   | 697 |    40 RETURN                                                             | 
                  
                          |   | 698 |       END                                                                | 
                  
                          |   | 699 |  | 
                  
                          |   | 700 | }}} | 
                  
                          |   | 701 |  | 
                  
                          |   | 702 | == Console19/Console19/Console19.f90 == | 
                  
                          |   | 703 | {{{ | 
                  
                          |   | 704 | #!fortran | 
                  
                          |   | 705 | !  Osc_rk_0.f90  | 
                  
                          |   | 706 | ! | 
                  
                          |   | 707 | !  Lab #3 | 
                  
                          |   | 708 | !  metod Runge-Kuta. Oscillator. | 
                  
                          |   | 709 |  | 
                  
                          |   | 710 | external fct, out | 
                  
                          |   | 711 | real aux(8,2) | 
                  
                          |   | 712 | common /a/ b, lam, f0, dt, k, n, c  | 
                  
                          |   | 713 |  | 
                  
                          |   | 714 | real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/) | 
                  
                          |   | 715 | real :: u(2) = (/0., 0.5/)  | 
                  
                          |   | 716 | real :: du(2) = (/0.5, 0.5/) | 
                  
                          |   | 717 | real :: b=0.003 | 
                  
                          |   | 718 | integer :: nd = 2 | 
                  
                          |   | 719 | real::  c = 2.0 | 
                  
                          |   | 720 | real :: dt = 2*3.14159265  | 
                  
                          |   | 721 |  | 
                  
                          |   | 722 |     | 
                  
                          |   | 723 | ! real :: lam = 0.1 | 
                  
                          |   | 724 | real :: f0 = 0.00 | 
                  
                          |   | 725 |  | 
                  
                          |   | 726 | n = 1  | 
                  
                          |   | 727 | k = 1 | 
                  
                          |   | 728 |   | 
                  
                          |   | 729 | open (unit=2, file = 'osc_rk.dat', status='new')  | 
                  
                          |   | 730 |  | 
                  
                          |   | 731 | write(*,*) 'x=', u(1), 'v=', u(2) | 
                  
                          |   | 732 | pause | 
                  
                          |   | 733 | call rkgs(pr, u, du, nd, ih, fct, out, aux) | 
                  
                          |   | 734 |  | 
                  
                          |   | 735 | write(2,*) 'ih=', ih | 
                  
                          |   | 736 |  | 
                  
                          |   | 737 | stop | 
                  
                          |   | 738 | end | 
                  
                          |   | 739 |  | 
                  
                          |   | 740 | subroutine fct (t, u, du) | 
                  
                          |   | 741 |  | 
                  
                          |   | 742 | real u(2), du(2), lam  | 
                  
                          |   | 743 | common /a/ b, lam, f0, dt, k, n, c | 
                  
                          |   | 744 |  | 
                  
                          |   | 745 | du(1) = c*cos(0.5*t)-u(2)-2.*b*u(1)                   | 
                  
                          |   | 746 | du(2) = u(1) | 
                  
                          |   | 747 |                                  | 
                  
                          |   | 748 | return  | 
                  
                          |   | 749 | end   | 
                  
                          |   | 750 |  | 
                  
                          |   | 751 | subroutine  out (t, u, du, ih, nd, pr) | 
                  
                          |   | 752 | real u(2), du(2), pr(5), lam | 
                  
                          |   | 753 | common /a/ b, lam, f0, dt, k, n, c | 
                  
                          |   | 754 |  | 
                  
                          |   | 755 | if (t.ge.k*(dt/20.)) then | 
                  
                          |   | 756 |  | 
                  
                          |   | 757 |  write (2,1) t,u(1),u(2)  !  | 
                  
                          |   | 758 |   k = k + 1 | 
                  
                          |   | 759 |    else | 
                  
                          |   | 760 |     end if | 
                  
                          |   | 761 |  | 
                  
                          |   | 762 | 1 format (f8.3, f8.3, f8.3) | 
                  
                          |   | 763 |  | 
                  
                          |   | 764 | return | 
                  
                          |   | 765 | end  | 
                  
                          |   | 766 |  | 
                  
                          |   | 767 |  | 
                  
                          |   | 768 |       SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)                | 
                  
                          |   | 769 | !                                                                        | 
                  
                          |   | 770 | !                                                                        | 
                  
                          |   | 771 |       DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)             | 
                  
                          |   | 772 |       DO 1 I=1,NDIM                                                      | 
                  
                          |   | 773 |     1 AUX(8,I)=.06666667*DERY(I)                                         | 
                  
                          |   | 774 |       X=PRMT(1)                                                          | 
                  
                          |   | 775 |       XEND=PRMT(2)                                                       | 
                  
                          |   | 776 |       H=PRMT(3)                                                          | 
                  
                          |   | 777 |       PRMT(5)=0.                                                         | 
                  
                          |   | 778 |       CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 779 | !                                                                        | 
                  
                          |   | 780 | !     ERROR TEST                                                         | 
                  
                          |   | 781 |       IF(H*(XEND-X))38,37,2                                              | 
                  
                          |   | 782 | !                                                                        | 
                  
                          |   | 783 | !     PREPARATIONS FOR RUNGE-KUTTA METHOD                                | 
                  
                          |   | 784 |     2 A(1)=.5                                                            | 
                  
                          |   | 785 |       A(2)=.2928932                                                      | 
                  
                          |   | 786 |       A(3)=1.707107                                                      | 
                  
                          |   | 787 |       A(4)=.1666667                                                      | 
                  
                          |   | 788 |       B(1)=2.                                                            | 
                  
                          |   | 789 |       B(2)=1.                                                            | 
                  
                          |   | 790 |       B(3)=1.                                                            | 
                  
                          |   | 791 |       B(4)=2.                                                            | 
                  
                          |   | 792 |       C(1)=.5                                                            | 
                  
                          |   | 793 |       C(2)=.2928932                                                      | 
                  
                          |   | 794 |       C(3)=1.707107                                                      | 
                  
                          |   | 795 |       C(4)=.5                                                            | 
                  
                          |   | 796 | !                                                                        | 
                  
                          |   | 797 | !     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                             | 
                  
                          |   | 798 |       DO 3 I=1,NDIM                                                      | 
                  
                          |   | 799 |       AUX(1,I)=Y(I)                                                      | 
                  
                          |   | 800 |       AUX(2,I)=DERY(I)                                                   | 
                  
                          |   | 801 |       AUX(3,I)=0.                                                        | 
                  
                          |   | 802 |     3 AUX(6,I)=0.                                                        | 
                  
                          |   | 803 |       IREC=0                                                             | 
                  
                          |   | 804 |       H=H+H                                                              | 
                  
                          |   | 805 |       IHLF=-1                                                            | 
                  
                          |   | 806 |       ISTEP=0                                                            | 
                  
                          |   | 807 |       IEND=0                                                             | 
                  
                          |   | 808 | !                                                                       | 
                  
                          |   | 809 | !                                                                    | 
                  
                          |   | 810 | !     START OF A RUNGE-KUTTA STEP                                        | 
                  
                          |   | 811 |     4 IF((X+H-XEND)*H)7,6,5                                              | 
                  
                          |   | 812 |     5 H=XEND-X                                                           | 
                  
                          |   | 813 |     6 IEND=1                                                             | 
                  
                          |   | 814 | !                                                                       | 
                  
                          |   | 815 | !     RECORDING OF INITIAL VALUES OF THIS STEP                           | 
                  
                          |   | 816 |     7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                 | 
                  
                          |   | 817 |       IF(PRMT(5))40,8,40                                                 | 
                  
                          |   | 818 |     8 ITEST=0                                                            | 
                  
                          |   | 819 |     9 ISTEP=ISTEP+1                                                      | 
                  
                          |   | 820 | !     START OF INNERMOST RUNGE-KUTTA LOOP                                | 
                  
                          |   | 821 |       J=1                                                                | 
                  
                          |   | 822 |    10 AJ=A(J)                                                            | 
                  
                          |   | 823 |       BJ=B(J)                                                            | 
                  
                          |   | 824 |       CJ=C(J)                                                            | 
                  
                          |   | 825 |       DO 11 I=1,NDIM                                                     | 
                  
                          |   | 826 |       R1=H*DERY(I)                                                       | 
                  
                          |   | 827 |       R2=AJ*(R1-BJ*AUX(6,I))                                             | 
                  
                          |   | 828 |       Y(I)=Y(I)+R2                                                       | 
                  
                          |   | 829 |       R2=R2+R2+R2                                                        | 
                  
                          |   | 830 |    11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                         | 
                  
                          |   | 831 |       IF(J-4)12,15,15                                                    | 
                  
                          |   | 832 |    12 J=J+1                                                              | 
                  
                          |   | 833 |       IF(J-3)13,14,13                                                    | 
                  
                          |   | 834 |    13 X=X+.5*H                                                           | 
                  
                          |   | 835 |    14 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 836 |       GOTO 10                                                            | 
                  
                          |   | 837 | !     END OF INNERMOST RUNGE-KUTTA LOOP                                  | 
                  
                          |   | 838 | !     TEST OF ACCURACY                                                   | 
                  
                          |   | 839 |    15 IF(ITEST)16,16,20                                                  | 
                  
                          |   | 840 | !                                                                        | 
                  
                          |   | 841 | !     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY    | 
                  
                          |   | 842 |    16 DO 17 I=1,NDIM                                                     | 
                  
                          |   | 843 |    17 AUX(4,I)=Y(I)                                                      | 
                  
                          |   | 844 |       ITEST=1                                                            | 
                  
                          |   | 845 |       ISTEP=ISTEP+ISTEP-2                                                | 
                  
                          |   | 846 |    18 IHLF=IHLF+1 | 
                  
                          |   | 847 |       X=X-H                                                              | 
                  
                          |   | 848 |       H=.5*H                                                             | 
                  
                          |   | 849 |       DO 19 I=1,NDIM                                                     | 
                  
                          |   | 850 |       Y(I)=AUX(1,I)                                                      | 
                  
                          |   | 851 |       DERY(I)=AUX(2,I)                                                   | 
                  
                          |   | 852 |    19 AUX(6,I)=AUX(3,I)                                                  | 
                  
                          |   | 853 |       GOTO 9                                                             | 
                  
                          |   | 854 | !                                                                        | 
                  
                          |   | 855 | !     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                    | 
                  
                          |   | 856 |    20 IMOD=ISTEP/2                                                       | 
                  
                          |   | 857 |       IF(ISTEP-IMOD-IMOD)21,23,21                                        | 
                  
                          |   | 858 |    21 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 859 |       DO 22 I=1,NDIM                                                     | 
                  
                          |   | 860 |       AUX(5,I)=Y(I)                                                      | 
                  
                          |   | 861 |    22 AUX(7,I)=DERY(I)                                                   | 
                  
                          |   | 862 |       GOTO 9                                                             | 
                  
                          |   | 863 | !                                                                        | 
                  
                          |   | 864 | !     COMPUTATION OF TEST VALUE DELT                                     | 
                  
                          |   | 865 |    23 DELT=0.                                                            | 
                  
                          |   | 866 |       DO 24 I=1,NDIM                                                     | 
                  
                          |   | 867 |    24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                              | 
                  
                          |   | 868 |       IF(DELT-PRMT(4))28,28,25                                           | 
                  
                          |   | 869 | !                                                                        | 
                  
                          |   | 870 | !     ERROR IS TOO GREAT                                                 | 
                  
                          |   | 871 |    25 IF(IHLF-10)26,36,36                                                | 
                  
                          |   | 872 |    26 DO 27 I=1,NDIM                                                     | 
                  
                          |   | 873 |    27 AUX(4,I)=AUX(5,I)                                                  | 
                  
                          |   | 874 |       ISTEP=ISTEP+ISTEP-4                                                | 
                  
                          |   | 875 |       X=X-H                                                              | 
                  
                          |   | 876 |       IEND=0                                                             | 
                  
                          |   | 877 |       GOTO 18                                                            | 
                  
                          |   | 878 | !                                                                        | 
                  
                          |   | 879 | !     RESULT VALUES ARE GOOD                                             | 
                  
                          |   | 880 |    28 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 881 |       DO 29 I=1,NDIM                                                     | 
                  
                          |   | 882 |       AUX(1,I)=Y(I)                                                      | 
                  
                          |   | 883 |       AUX(2,I)=DERY(I)                                                   | 
                  
                          |   | 884 |       AUX(3,I)=AUX(6,I)                                                  | 
                  
                          |   | 885 |       Y(I)=AUX(5,I)                                                      | 
                  
                          |   | 886 |    29 DERY(I)=AUX(7,I)                                                   | 
                  
                          |   | 887 |       CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                               | 
                  
                          |   | 888 |       IF(PRMT(5))40,30,40                                                | 
                  
                          |   | 889 |    30 DO 31 I=1,NDIM                                                     | 
                  
                          |   | 890 |       Y(I)=AUX(1,I)                                                      | 
                  
                          |   | 891 |    31 DERY(I)=AUX(2,I)                                                   | 
                  
                          |   | 892 |       IREC=IHLF                                                          | 
                  
                          |   | 893 |       IF(IEND)32,32,39                                                   | 
                  
                          |   | 894 | !                                                                        | 
                  
                          |   | 895 | !     INCREMENT GETS DOUBLED                                             | 
                  
                          |   | 896 |    32 IHLF=IHLF-1                                                        | 
                  
                          |   | 897 |       ISTEP=ISTEP/2                                                      | 
                  
                          |   | 898 |       H=H+H                                                              | 
                  
                          |   | 899 |       IF(IHLF)4,33,33                                                    | 
                  
                          |   | 900 |    33 IMOD=ISTEP/2                                                       | 
                  
                          |   | 901 |       IF(ISTEP-IMOD-IMOD)4,34,4                                          | 
                  
                          |   | 902 |    34 IF(DELT-.02*PRMT(4))35,35,4                                        | 
                  
                          |   | 903 |    35 IHLF=IHLF-1                                                        | 
                  
                          |   | 904 |       ISTEP=ISTEP/2                                                      | 
                  
                          |   | 905 |       H=H+H                                                              | 
                  
                          |   | 906 |       GOTO 4                                                             | 
                  
                          |   | 907 | !                                                                        | 
                  
                          |   | 908 | !     RETURNS TO CALLING PROGRAM                                         | 
                  
                          |   | 909 |    36 IHLF=11                                                            | 
                  
                          |   | 910 |       CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 911 |       GOTO 39                                                            | 
                  
                          |   | 912 |    37 IHLF=12                                                            | 
                  
                          |   | 913 |       GOTO 39                                                            | 
                  
                          |   | 914 |    38 IHLF=13                                                            | 
                  
                          |   | 915 |    39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                 | 
                  
                          |   | 916 |    40 RETURN                                                             | 
                  
                          |   | 917 |       END                                                                | 
                  
                          |   | 918 |  | 
                  
                          |   | 919 | }}} | 
                  
                          |   | 920 |  | 
                  
                          |   | 921 | == Console2/Console2/Console2.f90 == | 
                  
                          |   | 922 | {{{ | 
                  
                          |   | 923 | #!fortran | 
                  
                          |   | 924 | real x,dx,s | 
                  
                          |   | 925 | !open(1,file='y1.dat', status='unknown') | 
                  
                          |   | 926 | a=0 | 
                  
                          |   | 927 | b=1.57 | 
                  
                          |   | 928 | n=100 | 
                  
                          |   | 929 | dx=(b-a)/n | 
                  
                          |   | 930 | S=0 | 
                  
                          |   | 931 | do i=0,n | 
                  
                          |   | 932 |     if(mod(i,2)==0) then  | 
                  
                          |   | 933 |        S=S+dx/3*2*f(i*dx+a) | 
                  
                          |   | 934 |    else  | 
                  
                          |   | 935 |        S=S+dx/3*4*f(i*dx+a) | 
                  
                          |   | 936 |    end if | 
                  
                          |   | 937 | end do | 
                  
                          |   | 938 |   | 
                  
                          |   | 939 | print *, 'Sim = ', S  | 
                  
                          |   | 940 | pause | 
                  
                          |   | 941 |  | 
                  
                          |   | 942 | contains | 
                  
                          |   | 943 |  | 
                  
                          |   | 944 | function f(x)  | 
                  
                          |   | 945 | f=sin(x)*sin(2*x)*sin(3*x) | 
                  
                          |   | 946 | end function f | 
                  
                          |   | 947 | end | 
                  
                          |   | 948 | }}} | 
                  
                          |   | 949 |  | 
                  
                          |   | 950 | == Console20/Console20/Console20.f90 == | 
                  
                          |   | 951 | {{{ | 
                  
                          |   | 952 | #!fortran | 
                  
                          |   | 953 | !  Console20.f90  | 
                  
                          |   | 954 | ! | 
                  
                          |   | 955 | !  FUNCTIONS: | 
                  
                          |   | 956 | !  Console20 - Entry point of console application. | 
                  
                          |   | 957 | ! | 
                  
                          |   | 958 |  | 
                  
                          |   | 959 | !**************************************************************************** | 
                  
                          |   | 960 | ! | 
                  
                          |   | 961 | !  PROGRAM: Console20 | 
                  
                          |   | 962 | ! | 
                  
                          |   | 963 | !  PURPOSE:  Entry point for the console application. | 
                  
                          |   | 964 | ! | 
                  
                          |   | 965 | !**************************************************************************** | 
                  
                          |   | 966 |  | 
                  
                          |   | 967 |     program Console20 | 
                  
                          |   | 968 |  | 
                  
                          |   | 969 |     implicit none | 
                  
                          |   | 970 |  | 
                  
                          |   | 971 |     ! Variables | 
                  
                          |   | 972 |  | 
                  
                          |   | 973 |     ! Body of Console20 | 
                  
                          |   | 974 |     print *, 'Hello World' | 
                  
                          |   | 975 |  | 
                  
                          |   | 976 |     end program Console20 | 
                  
                          |   | 977 |  | 
                  
                          |   | 978 | }}} | 
                  
                          |   | 979 |  | 
                  
                          |   | 980 | == Console21/Console21/Console21.f90 == | 
                  
                          |   | 981 | {{{ | 
                  
                          |   | 982 | #!fortran | 
                  
                          |   | 983 | integer :: n=1000 | 
                  
                          |   | 984 | real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b | 
                  
                          |   | 985 | real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0  | 
                  
                          |   | 986 | x = 1.0 | 
                  
                          |   | 987 | y = 1.0 | 
                  
                          |   | 988 | a = 0.5 | 
                  
                          |   | 989 |  | 
                  
                          |   | 990 | vx = 1.0e09 | 
                  
                          |   | 991 | vy = 1.0e09 | 
                  
                          |   | 992 | E0 = 0.0 | 
                  
                          |   | 993 | B0 = 200.0 | 
                  
                          |   | 994 |  | 
                  
                          |   | 995 | b = (2*3.14)/(1.0e09) | 
                  
                          |   | 996 |  | 
                  
                          |   | 997 | t = 1.00e-08 | 
                  
                          |   | 998 | dt = t/float(n) | 
                  
                          |   | 999 |  | 
                  
                          |   | 1000 |  | 
                  
                          |   | 1001 | open (1,file='1',status='unknown') | 
                  
                          |   | 1002 |  | 
                  
                          |   | 1003 | do i=1,n | 
                  
                          |   | 1004 |     | 
                  
                          |   | 1005 |    Bz = B0*(1+a*y)  | 
                  
                          |   | 1006 |    !Bz = B0 | 
                  
                          |   | 1007 |    x = x + vx*dt | 
                  
                          |   | 1008 |    y = y + vy*dt | 
                  
                          |   | 1009 |  | 
                  
                          |   | 1010 |    vx = vx + (vy*eq*Bz)/(c*me)*dt | 
                  
                          |   | 1011 |    vy = vy - (vx*eq*Bz)/(c*me)*dt | 
                  
                          |   | 1012 |  | 
                  
                          |   | 1013 |     | 
                  
                          |   | 1014 |    if(i*dt.gt.b) then  | 
                  
                          |   | 1015 | E0 =  3.0 | 
                  
                          |   | 1016 |    end if | 
                  
                          |   | 1017 |     | 
                  
                          |   | 1018 |   vx = vx + (vy*eq*Bz)/(c*me)*dt | 
                  
                          |   | 1019 |   vy = vy + (eq*E0)/(me)*dt  - (vx*eq*Bz)/(c*me)*dt | 
                  
                          |   | 1020 |   x = x - vx*dt | 
                  
                          |   | 1021 |   y = y + vy*dt | 
                  
                          |   | 1022 |  | 
                  
                          |   | 1023 |    !print *, vx, vy | 
                  
                          |   | 1024 |  | 
                  
                          |   | 1025 |  | 
                  
                          |   | 1026 |  | 
                  
                          |   | 1027 |    | 
                  
                          |   | 1028 |  | 
                  
                          |   | 1029 |           write (1,*) x, y  | 
                  
                          |   | 1030 |  | 
                  
                          |   | 1031 |     end do | 
                  
                          |   | 1032 |  pause   | 
                  
                          |   | 1033 | end program | 
                  
                          |   | 1034 | }}} | 
                  
                          |   | 1035 |  | 
                  
                          |   | 1036 | == Console22/Console22/Console22.f90 == | 
                  
                          |   | 1037 | {{{ | 
                  
                          |   | 1038 | #!fortran | 
                  
                          |   | 1039 | integer :: n=1000 | 
                  
                          |   | 1040 | real(8) :: dt, t, x, y, vx, vy, v, E, W | 
                  
                          |   | 1041 | real(8) ::  mp = 1.67e-24, q = 4.8e-10 | 
                  
                          |   | 1042 | x = 1.0 | 
                  
                          |   | 1043 | y = 1.0 | 
                  
                          |   | 1044 | d = 1.0 | 
                  
                          |   | 1045 | a = 30 | 
                  
                          |   | 1046 | W0 = 1000 | 
                  
                          |   | 1047 |  | 
                  
                          |   | 1048 | v = sqrt(2*W0/mp) | 
                  
                          |   | 1049 |  | 
                  
                          |   | 1050 |  | 
                  
                          |   | 1051 | t = 1.00e-05 | 
                  
                          |   | 1052 | dt = t/float(n) | 
                  
                          |   | 1053 |  | 
                  
                          |   | 1054 |  | 
                  
                          |   | 1055 | open (1,file='1',status='unknown') | 
                  
                          |   | 1056 |  | 
                  
                          |   | 1057 | do i=1,n | 
                  
                          |   | 1058 |     | 
                  
                          |   | 1059 |   vx = v*sin(a) | 
                  
                          |   | 1060 |   vy = v*cos(a) | 
                  
                          |   | 1061 |   x = x + vx*dt | 
                  
                          |   | 1062 |    y = y + vy*dt | 
                  
                          |   | 1063 |    | 
                  
                          |   | 1064 |   E = (4*3.14*mp*y)/d | 
                  
                          |   | 1065 |   W = (E*q)/ x | 
                  
                          |   | 1066 |     | 
                  
                          |   | 1067 |  | 
                  
                          |   | 1068 |     | 
                  
                          |   | 1069 |  | 
                  
                          |   | 1070 |    !print *, vx, vy | 
                  
                          |   | 1071 |  | 
                  
                          |   | 1072 |  | 
                  
                          |   | 1073 |  | 
                  
                          |   | 1074 |    | 
                  
                          |   | 1075 |           write (1,*) W, t | 
                  
                          |   | 1076 |  | 
                  
                          |   | 1077 |     end do | 
                  
                          |   | 1078 |  pause   | 
                  
                          |   | 1079 | end program | 
                  
                          |   | 1080 | }}} | 
                  
                          |   | 1081 |  | 
                  
                          |   | 1082 | == Console24/Console24/Console24.f90 == | 
                  
                          |   | 1083 | {{{ | 
                  
                          |   | 1084 | #!fortran | 
                  
                          |   | 1085 |  | 
                  
                          |   | 1086 | !  Lab #1 | 
                  
                          |   | 1087 | !  metod Runge-Kuta. Rezonans i avtorezonans. | 
                  
                          |   | 1088 |  | 
                  
                          |   | 1089 | external fct, out | 
                  
                          |   | 1090 | real aux(8,2) | 
                  
                          |   | 1091 | common /a/  dt, k, n, g0, om, B0, B, a  | 
                  
                          |   | 1092 |  | 
                  
                          |   | 1093 | real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/) | 
                  
                          |   | 1094 | real :: u(2) = (/0., 0.5/)  | 
                  
                          |   | 1095 | real :: du(2) = (/0.5, 0.5/) | 
                  
                          |   | 1096 |  | 
                  
                          |   | 1097 | integer :: nd = 2 | 
                  
                          |   | 1098 | real :: dt = 2*3.14159265  | 
                  
                          |   | 1099 |  | 
                  
                          |   | 1100 | real :: t, m0, B, B0, a | 
                  
                          |   | 1101 | real :: f0 = 0.00 | 
                  
                          |   | 1102 | f = 2.45e09 | 
                  
                          |   | 1103 | pi = 3.14159265 | 
                  
                          |   | 1104 | ez = 4.8e-10 | 
                  
                          |   | 1105 | m0 = 9.1e-28 | 
                  
                          |   | 1106 | c = 3.0e10 | 
                  
                          |   | 1107 | om = 2*pi*f | 
                  
                          |   | 1108 | E = 2. | 
                  
                          |   | 1109 | w = 1000. | 
                  
                          |   | 1110 | a = 0.09 | 
                  
                          |   | 1111 | g0 = ez*E/(m0*c*om) | 
                  
                          |   | 1112 | B0 = m0*c*om/ez | 
                  
                          |   | 1113 | B = B0*(1+a*t) | 
                  
                          |   | 1114 | b = B/B0 | 
                  
                          |   | 1115 | u(1) = pi | 
                  
                          |   | 1116 | u(2) = w/511000.+1 | 
                  
                          |   | 1117 | write (*,*) g0 | 
                  
                          |   | 1118 | pause | 
                  
                          |   | 1119 |  | 
                  
                          |   | 1120 |  | 
                  
                          |   | 1121 |   | 
                  
                          |   | 1122 | n = 1  | 
                  
                          |   | 1123 | k = 1  | 
                  
                          |   | 1124 |   | 
                  
                          |   | 1125 | open (unit=2, file = 'rez_avtorez.dat', status='unknown')  | 
                  
                          |   | 1126 |  | 
                  
                          |   | 1127 |  | 
                  
                          |   | 1128 | call rkgs(pr, u, du, nd, ih, fct, out, aux) | 
                  
                          |   | 1129 |  | 
                  
                          |   | 1130 | write(2,*) 'ih=' ,ih  | 
                  
                          |   | 1131 |  | 
                  
                          |   | 1132 |  | 
                  
                          |   | 1133 |  | 
                  
                          |   | 1134 | end | 
                  
                          |   | 1135 |  | 
                  
                          |   | 1136 | subroutine fct (t, u, du) | 
                  
                          |   | 1137 |  | 
                  
                          |   | 1138 | real u(2), du(2), lam  | 
                  
                          |   | 1139 | common /a/ dt, k, n, g0, om, B, B0, a | 
                  
                          |   | 1140 |  | 
                  
                          |   | 1141 | du(1) = (a*t+1-u(2))/u(2)+g0*sqrt(1./(u(2)**2-1))*sin(u(1))           | 
                  
                          |   | 1142 | du(2) = -g0*sqrt(1.-1./u(2)**2)*cos(u(1)) | 
                  
                          |   | 1143 |  | 
                  
                          |   | 1144 | !write (*,*)  du(1), du(2) | 
                  
                          |   | 1145 | !pause | 
                  
                          |   | 1146 |                                  | 
                  
                          |   | 1147 | return  | 
                  
                          |   | 1148 | end   | 
                  
                          |   | 1149 |  | 
                  
                          |   | 1150 | subroutine  out (t, u, du, ih, nd, pr) | 
                  
                          |   | 1151 | real u(2), du(2), pr(5), lam | 
                  
                          |   | 1152 | common /a/ dt, k, n, g0, om, B, B0, a | 
                  
                          |   | 1153 |  | 
                  
                          |   | 1154 | if (t.ge.k*(dt/2.)) then | 
                  
                          |   | 1155 |  | 
                  
                          |   | 1156 |  write (2,1) t,u(1),(u(2)-1)*511.  | 
                  
                          |   | 1157 |   k = k + 1 | 
                  
                          |   | 1158 |    else | 
                  
                          |   | 1159 |     end if | 
                  
                          |   | 1160 |  | 
                  
                          |   | 1161 | 1 format (3e10.3) | 
                  
                          |   | 1162 |  | 
                  
                          |   | 1163 | return | 
                  
                          |   | 1164 | end  | 
                  
                          |   | 1165 |  | 
                  
                          |   | 1166 |  | 
                  
                          |   | 1167 |       SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)                | 
                  
                          |   | 1168 | !                                                                        | 
                  
                          |   | 1169 | !                                                                        | 
                  
                          |   | 1170 |       DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)             | 
                  
                          |   | 1171 |       DO 1 I=1,NDIM                                                      | 
                  
                          |   | 1172 |     1 AUX(8,I)=.06666667*DERY(I)                                         | 
                  
                          |   | 1173 |       X=PRMT(1)                                                          | 
                  
                          |   | 1174 |       XEND=PRMT(2)                                                       | 
                  
                          |   | 1175 |       H=PRMT(3)                                                          | 
                  
                          |   | 1176 |       PRMT(5)=0.                                                         | 
                  
                          |   | 1177 |       CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 1178 | !                                                                        | 
                  
                          |   | 1179 | !     ERROR TEST                                                         | 
                  
                          |   | 1180 |       IF(H*(XEND-X))38,37,2                                              | 
                  
                          |   | 1181 | !                                                                        | 
                  
                          |   | 1182 | !     PREPARATIONS FOR RUNGE-KUTTA METHOD                                | 
                  
                          |   | 1183 |     2 A(1)=.5                                                            | 
                  
                          |   | 1184 |       A(2)=.2928932                                                      | 
                  
                          |   | 1185 |       A(3)=1.707107                                                      | 
                  
                          |   | 1186 |       A(4)=.1666667                                                      | 
                  
                          |   | 1187 |       B(1)=2.                                                            | 
                  
                          |   | 1188 |       B(2)=1.                                                            | 
                  
                          |   | 1189 |       B(3)=1.                                                            | 
                  
                          |   | 1190 |       B(4)=2.                                                            | 
                  
                          |   | 1191 |       C(1)=.5                                                            | 
                  
                          |   | 1192 |       C(2)=.2928932                                                      | 
                  
                          |   | 1193 |       C(3)=1.707107                                                      | 
                  
                          |   | 1194 |       C(4)=.5                                                            | 
                  
                          |   | 1195 | !                                                                        | 
                  
                          |   | 1196 | !     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                             | 
                  
                          |   | 1197 |       DO 3 I=1,NDIM                                                      | 
                  
                          |   | 1198 |       AUX(1,I)=Y(I)                                                      | 
                  
                          |   | 1199 |       AUX(2,I)=DERY(I)                                                   | 
                  
                          |   | 1200 |       AUX(3,I)=0.                                                        | 
                  
                          |   | 1201 |     3 AUX(6,I)=0.                                                        | 
                  
                          |   | 1202 |       IREC=0                                                             | 
                  
                          |   | 1203 |       H=H+H                                                              | 
                  
                          |   | 1204 |       IHLF=-1                                                            | 
                  
                          |   | 1205 |       ISTEP=0                                                            | 
                  
                          |   | 1206 |       IEND=0                                                             | 
                  
                          |   | 1207 | !                                                                       | 
                  
                          |   | 1208 | !                                                                    | 
                  
                          |   | 1209 | !     START OF A RUNGE-KUTTA STEP                                        | 
                  
                          |   | 1210 |     4 IF((X+H-XEND)*H)7,6,5                                              | 
                  
                          |   | 1211 |     5 H=XEND-X                                                           | 
                  
                          |   | 1212 |     6 IEND=1                                                             | 
                  
                          |   | 1213 | !                                                                       | 
                  
                          |   | 1214 | !     RECORDING OF INITIAL VALUES OF THIS STEP                           | 
                  
                          |   | 1215 |     7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                 | 
                  
                          |   | 1216 |       IF(PRMT(5))40,8,40                                                 | 
                  
                          |   | 1217 |     8 ITEST=0                                                            | 
                  
                          |   | 1218 |     9 ISTEP=ISTEP+1                                                      | 
                  
                          |   | 1219 | !     START OF INNERMOST RUNGE-KUTTA LOOP                                | 
                  
                          |   | 1220 |       J=1                                                                | 
                  
                          |   | 1221 |    10 AJ=A(J)                                                            | 
                  
                          |   | 1222 |       BJ=B(J)                                                            | 
                  
                          |   | 1223 |       CJ=C(J)                                                            | 
                  
                          |   | 1224 |       DO 11 I=1,NDIM                                                     | 
                  
                          |   | 1225 |       R1=H*DERY(I)                                                       | 
                  
                          |   | 1226 |       R2=AJ*(R1-BJ*AUX(6,I))                                             | 
                  
                          |   | 1227 |       Y(I)=Y(I)+R2                                                       | 
                  
                          |   | 1228 |       R2=R2+R2+R2                                                        | 
                  
                          |   | 1229 |    11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                         | 
                  
                          |   | 1230 |       IF(J-4)12,15,15                                                    | 
                  
                          |   | 1231 |    12 J=J+1                                                              | 
                  
                          |   | 1232 |       IF(J-3)13,14,13                                                    | 
                  
                          |   | 1233 |    13 X=X+.5*H                                                           | 
                  
                          |   | 1234 |    14 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 1235 |       GOTO 10                                                            | 
                  
                          |   | 1236 | !     END OF INNERMOST RUNGE-KUTTA LOOP                                  | 
                  
                          |   | 1237 | !     TEST OF ACCURACY                                                   | 
                  
                          |   | 1238 |    15 IF(ITEST)16,16,20                                                  | 
                  
                          |   | 1239 | !                                                                        | 
                  
                          |   | 1240 | !     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY    | 
                  
                          |   | 1241 |    16 DO 17 I=1,NDIM                                                     | 
                  
                          |   | 1242 |    17 AUX(4,I)=Y(I)                                                      | 
                  
                          |   | 1243 |       ITEST=1                                                            | 
                  
                          |   | 1244 |       ISTEP=ISTEP+ISTEP-2                                                | 
                  
                          |   | 1245 |    18 IHLF=IHLF+1 | 
                  
                          |   | 1246 |       X=X-H                                                              | 
                  
                          |   | 1247 |       H=.5*H                                                             | 
                  
                          |   | 1248 |       DO 19 I=1,NDIM                                                     | 
                  
                          |   | 1249 |       Y(I)=AUX(1,I)                                                      | 
                  
                          |   | 1250 |       DERY(I)=AUX(2,I)                                                   | 
                  
                          |   | 1251 |    19 AUX(6,I)=AUX(3,I)                                                  | 
                  
                          |   | 1252 |       GOTO 9                                                             | 
                  
                          |   | 1253 | !                                                                        | 
                  
                          |   | 1254 | !     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                    | 
                  
                          |   | 1255 |    20 IMOD=ISTEP/2                                                       | 
                  
                          |   | 1256 |       IF(ISTEP-IMOD-IMOD)21,23,21                                        | 
                  
                          |   | 1257 |    21 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 1258 |       DO 22 I=1,NDIM                                                     | 
                  
                          |   | 1259 |       AUX(5,I)=Y(I)                                                      | 
                  
                          |   | 1260 |    22 AUX(7,I)=DERY(I)                                                   | 
                  
                          |   | 1261 |       GOTO 9                                                             | 
                  
                          |   | 1262 | !                                                                        | 
                  
                          |   | 1263 | !     COMPUTATION OF TEST VALUE DELT                                     | 
                  
                          |   | 1264 |    23 DELT=0.                                                            | 
                  
                          |   | 1265 |       DO 24 I=1,NDIM                                                     | 
                  
                          |   | 1266 |    24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                              | 
                  
                          |   | 1267 |       IF(DELT-PRMT(4))28,28,25                                           | 
                  
                          |   | 1268 | !                                                                        | 
                  
                          |   | 1269 | !     ERROR IS TOO GREAT                                                 | 
                  
                          |   | 1270 |    25 IF(IHLF-10)26,36,36                                                | 
                  
                          |   | 1271 |    26 DO 27 I=1,NDIM                                                     | 
                  
                          |   | 1272 |    27 AUX(4,I)=AUX(5,I)                                                  | 
                  
                          |   | 1273 |       ISTEP=ISTEP+ISTEP-4                                                | 
                  
                          |   | 1274 |       X=X-H                                                              | 
                  
                          |   | 1275 |       IEND=0                                                             | 
                  
                          |   | 1276 |       GOTO 18                                                            | 
                  
                          |   | 1277 | !                                                                        | 
                  
                          |   | 1278 | !     RESULT VALUES ARE GOOD                                             | 
                  
                          |   | 1279 |    28 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 1280 |       DO 29 I=1,NDIM                                                     | 
                  
                          |   | 1281 |       AUX(1,I)=Y(I)                                                      | 
                  
                          |   | 1282 |       AUX(2,I)=DERY(I)                                                   | 
                  
                          |   | 1283 |       AUX(3,I)=AUX(6,I)                                                  | 
                  
                          |   | 1284 |       Y(I)=AUX(5,I)                                                      | 
                  
                          |   | 1285 |    29 DERY(I)=AUX(7,I)                                                   | 
                  
                          |   | 1286 |       CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                               | 
                  
                          |   | 1287 |       IF(PRMT(5))40,30,40                                                | 
                  
                          |   | 1288 |    30 DO 31 I=1,NDIM                                                     | 
                  
                          |   | 1289 |       Y(I)=AUX(1,I)                                                      | 
                  
                          |   | 1290 |    31 DERY(I)=AUX(2,I)                                                   | 
                  
                          |   | 1291 |       IREC=IHLF                                                          | 
                  
                          |   | 1292 |       IF(IEND)32,32,39                                                   | 
                  
                          |   | 1293 | !                                                                        | 
                  
                          |   | 1294 | !     INCREMENT GETS DOUBLED                                             | 
                  
                          |   | 1295 |    32 IHLF=IHLF-1                                                        | 
                  
                          |   | 1296 |       ISTEP=ISTEP/2                                                      | 
                  
                          |   | 1297 |       H=H+H                                                              | 
                  
                          |   | 1298 |       IF(IHLF)4,33,33                                                    | 
                  
                          |   | 1299 |    33 IMOD=ISTEP/2                                                       | 
                  
                          |   | 1300 |       IF(ISTEP-IMOD-IMOD)4,34,4                                          | 
                  
                          |   | 1301 |    34 IF(DELT-.02*PRMT(4))35,35,4                                        | 
                  
                          |   | 1302 |    35 IHLF=IHLF-1                                                        | 
                  
                          |   | 1303 |       ISTEP=ISTEP/2                                                      | 
                  
                          |   | 1304 |       H=H+H                                                              | 
                  
                          |   | 1305 |       GOTO 4                                                             | 
                  
                          |   | 1306 | !                                                                        | 
                  
                          |   | 1307 | !     RETURNS TO CALLING PROGRAM                                         | 
                  
                          |   | 1308 |    36 IHLF=11                                                            | 
                  
                          |   | 1309 |       CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 1310 |       GOTO 39                                                            | 
                  
                          |   | 1311 |    37 IHLF=12                                                            | 
                  
                          |   | 1312 |       GOTO 39                                                            | 
                  
                          |   | 1313 |    38 IHLF=13                                                            | 
                  
                          |   | 1314 |    39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                 | 
                  
                          |   | 1315 |    40 RETURN                                                             | 
                  
                          |   | 1316 |       END                                                                | 
                  
                          |   | 1317 |  | 
                  
                          |   | 1318 |  | 
                  
                          |   | 1319 | }}} | 
                  
                          |   | 1320 |  | 
                  
                          |   | 1321 | == Console26/Console26/Console26.f90 == | 
                  
                          |   | 1322 | {{{ | 
                  
                          |   | 1323 | #!fortran | 
                  
                          |   | 1324 |  | 
                  
                          |   | 1325 | !  Lab #1 | 
                  
                          |   | 1326 | !  metod Runge-Kuta. Rezonans i avtorezonans. | 
                  
                          |   | 1327 |  | 
                  
                          |   | 1328 | external fct, out | 
                  
                          |   | 1329 | real aux(8,4) | 
                  
                          |   | 1330 | common /a/  dt, k, n, g0, om, B0, B, a, y, rl  | 
                  
                          |   | 1331 |  | 
                  
                          |   | 1332 | real :: pr(5)=(/0., 2*3.14159*200., 2*3.14159/50, .0001, .0/) | 
                  
                          |   | 1333 | real :: u(4) = (/0., 0., 0., 0./)  | 
                  
                          |   | 1334 | real :: du(4) = (/0.25, 0.25, 0.25, 0.25/) | 
                  
                          |   | 1335 |  | 
                  
                          |   | 1336 | integer :: nd = 4 | 
                  
                          |   | 1337 | real :: dt = 2*3.14159265  | 
                  
                          |   | 1338 |  | 
                  
                          |   | 1339 | real :: t, m0, B, B0, a | 
                  
                          |   | 1340 | real :: f0 = 0.00 | 
                  
                          |   | 1341 | f = 2.45e09 | 
                  
                          |   | 1342 | pi = 3.14159265 | 
                  
                          |   | 1343 | ez = 4.8e-10 | 
                  
                          |   | 1344 | m0 = 9.1e-28 | 
                  
                          |   | 1345 | c = 3.0e10 | 
                  
                          |   | 1346 | om = 2*pi*f | 
                  
                          |   | 1347 | E0 = 5. | 
                  
                          |   | 1348 | E = E0 | 
                  
                          |   | 1349 | a = 0.00034 | 
                  
                          |   | 1350 | B0 = m0*c*om/ez | 
                  
                          |   | 1351 | g0 = E0/B0 | 
                  
                          |   | 1352 | B = B0 | 
                  
                          |   | 1353 | b = B/B0 | 
                  
                          |   | 1354 | rl = c/om | 
                  
                          |   | 1355 | u(1) = 0. | 
                  
                          |   | 1356 | u(2) = 0. | 
                  
                          |   | 1357 | u(3) = 0. | 
                  
                          |   | 1358 | u(4) = 0. | 
                  
                          |   | 1359 | write (*,*) g0, B0, b, E, rl | 
                  
                          |   | 1360 | pause | 
                  
                          |   | 1361 |  | 
                  
                          |   | 1362 |  | 
                  
                          |   | 1363 |   | 
                  
                          |   | 1364 | n = 1  | 
                  
                          |   | 1365 | k = 1  | 
                  
                          |   | 1366 |   | 
                  
                          |   | 1367 | open (unit=2, file = 'rez_avtorez.dat', status='unknown')  | 
                  
                          |   | 1368 |  | 
                  
                          |   | 1369 |  | 
                  
                          |   | 1370 | call rkgs(pr, u, du, nd, ih, fct, out, aux) | 
                  
                          |   | 1371 |  | 
                  
                          |   | 1372 | write(2,*) 'ih=' ,ih  | 
                  
                          |   | 1373 |  | 
                  
                          |   | 1374 |  | 
                  
                          |   | 1375 |  | 
                  
                          |   | 1376 | end | 
                  
                          |   | 1377 |  | 
                  
                          |   | 1378 | subroutine fct (t, u, du) | 
                  
                          |   | 1379 |  | 
                  
                          |   | 1380 | real u(4), du(4), lam  | 
                  
                          |   | 1381 | common /a/ dt, k, n, g0, om, B, B0, a, y, rl | 
                  
                          |   | 1382 | y=sqrt(1+u(1)**2+u(2)**2) | 
                  
                          |   | 1383 | du(1) = -g0*cos(t)- u(2)*(1+a*t)/y        | 
                  
                          |   | 1384 | du(2) = -g0*sin(t)+u(1)*(1+a*t)/y | 
                  
                          |   | 1385 | du(3) = u(1)/y | 
                  
                          |   | 1386 | du(4) = u(2)/y | 
                  
                          |   | 1387 |  | 
                  
                          |   | 1388 | !write (*,*)  du(1), du(2), du(3), du(4) | 
                  
                          |   | 1389 | !pause | 
                  
                          |   | 1390 |                                  | 
                  
                          |   | 1391 | return  | 
                  
                          |   | 1392 | end   | 
                  
                          |   | 1393 |  | 
                  
                          |   | 1394 | subroutine  out (t, u, du, ih, nd, pr) | 
                  
                          |   | 1395 | real u(4), du(4), pr(5), lam | 
                  
                          |   | 1396 | common /a/ dt, k, n, g0, om, B, B0, a, y, rl | 
                  
                          |   | 1397 |  | 
                  
                          |   | 1398 | if (t.ge.k*(dt)) then | 
                  
                          |   | 1399 |  | 
                  
                          |   | 1400 |  write (2,1) t,(y-1)*511. ,u(3)*rl,u(4)*rl  | 
                  
                          |   | 1401 |   k = k + 1 | 
                  
                          |   | 1402 |    else | 
                  
                          |   | 1403 |     end if | 
                  
                          |   | 1404 |  | 
                  
                          |   | 1405 | 1 format (5e12.3) | 
                  
                          |   | 1406 |  | 
                  
                          |   | 1407 | return | 
                  
                          |   | 1408 | end  | 
                  
                          |   | 1409 |  | 
                  
                          |   | 1410 |  | 
                  
                          |   | 1411 |       SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)                | 
                  
                          |   | 1412 | !                                                                        | 
                  
                          |   | 1413 | !                                                                        | 
                  
                          |   | 1414 |       DIMENSION Y(4),DERY(4),AUX(8,4),A(4),B(4),C(4),PRMT(5)             | 
                  
                          |   | 1415 |       DO 1 I=1,NDIM                                                      | 
                  
                          |   | 1416 |     1 AUX(8,I)=.06666667*DERY(I)                                         | 
                  
                          |   | 1417 |       X=PRMT(1)                                                          | 
                  
                          |   | 1418 |       XEND=PRMT(2)                                                       | 
                  
                          |   | 1419 |       H=PRMT(3)                                                          | 
                  
                          |   | 1420 |       PRMT(5)=0.                                                         | 
                  
                          |   | 1421 |       CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 1422 | !                                                                        | 
                  
                          |   | 1423 | !     ERROR TEST                                                         | 
                  
                          |   | 1424 |       IF(H*(XEND-X))38,37,2                                              | 
                  
                          |   | 1425 | !                                                                        | 
                  
                          |   | 1426 | !     PREPARATIONS FOR RUNGE-KUTTA METHOD                                | 
                  
                          |   | 1427 |     2 A(1)=.5                                                            | 
                  
                          |   | 1428 |       A(2)=.2928932                                                      | 
                  
                          |   | 1429 |       A(3)=1.707107                                                      | 
                  
                          |   | 1430 |       A(4)=.1666667                                                      | 
                  
                          |   | 1431 |       B(1)=2.                                                            | 
                  
                          |   | 1432 |       B(2)=1.                                                            | 
                  
                          |   | 1433 |       B(3)=1.                                                            | 
                  
                          |   | 1434 |       B(4)=2.                                                            | 
                  
                          |   | 1435 |       C(1)=.5                                                            | 
                  
                          |   | 1436 |       C(2)=.2928932                                                      | 
                  
                          |   | 1437 |       C(3)=1.707107                                                      | 
                  
                          |   | 1438 |       C(4)=.5                                                            | 
                  
                          |   | 1439 | !                                                                        | 
                  
                          |   | 1440 | !     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                             | 
                  
                          |   | 1441 |       DO 3 I=1,NDIM                                                      | 
                  
                          |   | 1442 |       AUX(1,I)=Y(I)                                                      | 
                  
                          |   | 1443 |       AUX(2,I)=DERY(I)                                                   | 
                  
                          |   | 1444 |       AUX(3,I)=0.                                                        | 
                  
                          |   | 1445 |     3 AUX(6,I)=0.                                                        | 
                  
                          |   | 1446 |       IREC=0                                                             | 
                  
                          |   | 1447 |       H=H+H                                                              | 
                  
                          |   | 1448 |       IHLF=-1                                                            | 
                  
                          |   | 1449 |       ISTEP=0                                                            | 
                  
                          |   | 1450 |       IEND=0                                                             | 
                  
                          |   | 1451 | !                                                                       | 
                  
                          |   | 1452 | !                                                                    | 
                  
                          |   | 1453 | !     START OF A RUNGE-KUTTA STEP                                        | 
                  
                          |   | 1454 |     4 IF((X+H-XEND)*H)7,6,5                                              | 
                  
                          |   | 1455 |     5 H=XEND-X                                                           | 
                  
                          |   | 1456 |     6 IEND=1                                                             | 
                  
                          |   | 1457 | !                                                                       | 
                  
                          |   | 1458 | !     RECORDING OF INITIAL VALUES OF THIS STEP                           | 
                  
                          |   | 1459 |     7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                 | 
                  
                          |   | 1460 |       IF(PRMT(5))40,8,40                                                 | 
                  
                          |   | 1461 |     8 ITEST=0                                                            | 
                  
                          |   | 1462 |     9 ISTEP=ISTEP+1                                                      | 
                  
                          |   | 1463 | !     START OF INNERMOST RUNGE-KUTTA LOOP                                | 
                  
                          |   | 1464 |       J=1                                                                | 
                  
                          |   | 1465 |    10 AJ=A(J)                                                            | 
                  
                          |   | 1466 |       BJ=B(J)                                                            | 
                  
                          |   | 1467 |       CJ=C(J)                                                            | 
                  
                          |   | 1468 |       DO 11 I=1,NDIM                                                     | 
                  
                          |   | 1469 |       R1=H*DERY(I)                                                       | 
                  
                          |   | 1470 |       R2=AJ*(R1-BJ*AUX(6,I))                                             | 
                  
                          |   | 1471 |       Y(I)=Y(I)+R2                                                       | 
                  
                          |   | 1472 |       R2=R2+R2+R2                                                        | 
                  
                          |   | 1473 |    11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                         | 
                  
                          |   | 1474 |       IF(J-4)12,15,15                                                    | 
                  
                          |   | 1475 |    12 J=J+1                                                              | 
                  
                          |   | 1476 |       IF(J-3)13,14,13                                                    | 
                  
                          |   | 1477 |    13 X=X+.5*H                                                           | 
                  
                          |   | 1478 |    14 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 1479 |       GOTO 10                                                            | 
                  
                          |   | 1480 | !     END OF INNERMOST RUNGE-KUTTA LOOP                                  | 
                  
                          |   | 1481 | !     TEST OF ACCURACY                                                   | 
                  
                          |   | 1482 |    15 IF(ITEST)16,16,20                                                  | 
                  
                          |   | 1483 | !                                                                        | 
                  
                          |   | 1484 | !     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY    | 
                  
                          |   | 1485 |    16 DO 17 I=1,NDIM                                                     | 
                  
                          |   | 1486 |    17 AUX(4,I)=Y(I)                                                      | 
                  
                          |   | 1487 |       ITEST=1                                                            | 
                  
                          |   | 1488 |       ISTEP=ISTEP+ISTEP-2                                                | 
                  
                          |   | 1489 |    18 IHLF=IHLF+1 | 
                  
                          |   | 1490 |       X=X-H                                                              | 
                  
                          |   | 1491 |       H=.5*H                                                             | 
                  
                          |   | 1492 |       DO 19 I=1,NDIM                                                     | 
                  
                          |   | 1493 |       Y(I)=AUX(1,I)                                                      | 
                  
                          |   | 1494 |       DERY(I)=AUX(2,I)                                                   | 
                  
                          |   | 1495 |    19 AUX(6,I)=AUX(3,I)                                                  | 
                  
                          |   | 1496 |       GOTO 9                                                             | 
                  
                          |   | 1497 | !                                                                        | 
                  
                          |   | 1498 | !     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                    | 
                  
                          |   | 1499 |    20 IMOD=ISTEP/2                                                       | 
                  
                          |   | 1500 |       IF(ISTEP-IMOD-IMOD)21,23,21                                        | 
                  
                          |   | 1501 |    21 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 1502 |       DO 22 I=1,NDIM                                                     | 
                  
                          |   | 1503 |       AUX(5,I)=Y(I)                                                      | 
                  
                          |   | 1504 |    22 AUX(7,I)=DERY(I)                                                   | 
                  
                          |   | 1505 |       GOTO 9                                                             | 
                  
                          |   | 1506 | !                                                                        | 
                  
                          |   | 1507 | !     COMPUTATION OF TEST VALUE DELT                                     | 
                  
                          |   | 1508 |    23 DELT=0.                                                            | 
                  
                          |   | 1509 |       DO 24 I=1,NDIM                                                     | 
                  
                          |   | 1510 |    24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                              | 
                  
                          |   | 1511 |       IF(DELT-PRMT(4))28,28,25                                           | 
                  
                          |   | 1512 | !                                                                        | 
                  
                          |   | 1513 | !     ERROR IS TOO GREAT                                                 | 
                  
                          |   | 1514 |    25 IF(IHLF-10)26,36,36                                                | 
                  
                          |   | 1515 |    26 DO 27 I=1,NDIM                                                     | 
                  
                          |   | 1516 |    27 AUX(4,I)=AUX(5,I)                                                  | 
                  
                          |   | 1517 |       ISTEP=ISTEP+ISTEP-4                                                | 
                  
                          |   | 1518 |       X=X-H                                                              | 
                  
                          |   | 1519 |       IEND=0                                                             | 
                  
                          |   | 1520 |       GOTO 18                                                            | 
                  
                          |   | 1521 | !                                                                        | 
                  
                          |   | 1522 | !     RESULT VALUES ARE GOOD                                             | 
                  
                          |   | 1523 |    28 CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 1524 |       DO 29 I=1,NDIM                                                     | 
                  
                          |   | 1525 |       AUX(1,I)=Y(I)                                                      | 
                  
                          |   | 1526 |       AUX(2,I)=DERY(I)                                                   | 
                  
                          |   | 1527 |       AUX(3,I)=AUX(6,I)                                                  | 
                  
                          |   | 1528 |       Y(I)=AUX(5,I)                                                      | 
                  
                          |   | 1529 |    29 DERY(I)=AUX(7,I)                                                   | 
                  
                          |   | 1530 |       CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                               | 
                  
                          |   | 1531 |       IF(PRMT(5))40,30,40                                                | 
                  
                          |   | 1532 |    30 DO 31 I=1,NDIM                                                     | 
                  
                          |   | 1533 |       Y(I)=AUX(1,I)                                                      | 
                  
                          |   | 1534 |    31 DERY(I)=AUX(2,I)                                                   | 
                  
                          |   | 1535 |       IREC=IHLF                                                          | 
                  
                          |   | 1536 |       IF(IEND)32,32,39                                                   | 
                  
                          |   | 1537 | !                                                                        | 
                  
                          |   | 1538 | !     INCREMENT GETS DOUBLED                                             | 
                  
                          |   | 1539 |    32 IHLF=IHLF-1                                                        | 
                  
                          |   | 1540 |       ISTEP=ISTEP/2                                                      | 
                  
                          |   | 1541 |       H=H+H                                                              | 
                  
                          |   | 1542 |       IF(IHLF)4,33,33                                                    | 
                  
                          |   | 1543 |    33 IMOD=ISTEP/2                                                       | 
                  
                          |   | 1544 |       IF(ISTEP-IMOD-IMOD)4,34,4                                          | 
                  
                          |   | 1545 |    34 IF(DELT-.02*PRMT(4))35,35,4                                        | 
                  
                          |   | 1546 |    35 IHLF=IHLF-1                                                        | 
                  
                          |   | 1547 |       ISTEP=ISTEP/2                                                      | 
                  
                          |   | 1548 |       H=H+H                                                              | 
                  
                          |   | 1549 |       GOTO 4                                                             | 
                  
                          |   | 1550 | !                                                                        | 
                  
                          |   | 1551 | !     RETURNS TO CALLING PROGRAM                                         | 
                  
                          |   | 1552 |    36 IHLF=11                                                            | 
                  
                          |   | 1553 |       CALL FCT(X,Y,DERY)                                                 | 
                  
                          |   | 1554 |       GOTO 39       | 
                  
                          |   | 1555 |    38 IHLF=13                                                            | 
                  
                          |   | 1556 |    39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                 | 
                  
                          |   | 1557 |    40 RETURN                                                             | 
                  
                          |   | 1558 |       END                                                                | 
                  
                          |   | 1559 | }}} | 
                  
                          |   | 1560 |    |