home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / SHOOTF.DEM < prev    next >
Encoding:
Text File  |  1991-04-29  |  3.3 KB  |  143 lines
  1. PROGRAM d16r2(input,output);
  2. (* driver for routine SHOOTF *)
  3. LABEL 1;
  4. CONST
  5.    nvar=3;
  6.    n1=2;
  7.    n2=1;
  8.    delta=1.0e-3;
  9.    eps=1.0e-6;
  10.    dx=1.0e-4;
  11. TYPE
  12.    gl1array = ARRAY [1..1] OF real;
  13.    gln2array = gl1array;
  14.    gl2array = ARRAY [1..2] OF real;
  15.    gln1array = gl2array;
  16.    gl3array = ARRAY [1..3] OF real;
  17.    glnarray = gl3array;
  18.    glarray = gl3array;
  19.    glnvar = gl3array;
  20.    glnpbynp = ARRAY [1..3,1..3] OF real;
  21.    glnvarbynvar = glnpbynp;
  22.    glindx = ARRAY [1..3] OF integer;
  23. VAR
  24.    c2,factr,h1,hmin : real;
  25.    q1,x1,x2,xf : real;
  26.    i,m,n : integer;
  27.    v1,delv1,dv1 : gln2array;
  28.    v2,delv2,dv2 : gln1array;
  29.    f : glnvar;
  30.    kmax,kount : integer;
  31.    dxsav : real;
  32.    xp : ARRAY [1..200] OF real;
  33.    yp : ARRAY [1..10,1..200] OF real;
  34.  
  35. PROCEDURE load1(x1: real; v1: gl1array; VAR y: gl3array);
  36. (* Programs using routine LOAD1 must define the types
  37. TYPE
  38.    gl1array = ARRAY [1..1] OF real;
  39.    gl3array = ARRAY [1..3] OF real;
  40. and must declare the variables
  41. VAR
  42.    c2,factr : real;
  43.    m : integer;
  44. in the main routine. *)
  45. BEGIN
  46.    y[3] := v1[1];
  47.    y[2] := -(y[3]-c2)*factr/2.0/(m+1.0);
  48.    y[1] := factr+y[2]*dx
  49. END;
  50.  
  51. PROCEDURE load2(x2: real; v2: gl2array; VAR y: gl3array);
  52. (* Programs using routine LOAD2 must define the types
  53. TYPE
  54.    gl2array = ARRAY [1..2] OF real;
  55.    gl3array = ARRAY [1..3] OF real;
  56. and must declare the variables
  57.    c2 : real;
  58.    m : integer;
  59. in the main routine. *)
  60. BEGIN
  61.    y[3] := v2[2];
  62.    y[1] := v2[1];
  63.    y[2] := (y[3]-c2)*y[1]/2.0/(m+1.0)
  64. END;
  65.  
  66. PROCEDURE score(xf: real; y: gl3array; VAR f: gl3array);
  67. (* Programs using routine SCORE must define the type
  68. TYPE
  69.    gl3array = ARRAY [1..3] OF real;
  70. in the main routine. *)
  71. VAR
  72.    i : integer;
  73. BEGIN
  74.    FOR i := 1 to 3 DO BEGIN
  75.       f[i] := y[i]
  76.    END
  77. END;
  78.  
  79. PROCEDURE derivs(x: real; y: gl3array; VAR dydx: gl3array);
  80. (* Programs using routine DERIVS must define the type
  81. TYPE
  82.    gl3array = ARRAY [1..3] OF real;
  83. and must declare the variables
  84.    c2 : real;
  85.    m : integer;
  86. in the main routine. *)
  87. BEGIN
  88.    dydx[1] := y[2];
  89.    dydx[3] := 0.0;
  90.    dydx[2] := (2.0*x*(m+1.0)*y[2]-(y[3]-c2*x*x)*y[1])/(1.0-x*x)
  91. END;
  92.  
  93. (*$I MODFILE.PAS *)
  94. (*$I LUBKSB.PAS *)
  95.  
  96. (*$I LUDCMP.PAS *)
  97.  
  98. (*$I RK4.PAS *)
  99.  
  100. (*$I RKQC.PAS *)
  101.  
  102. (*$I ODEINT.PAS *)
  103.  
  104. (*$I SHOOTF.PAS *)
  105.  
  106. BEGIN
  107. 1:   write('Input M,N,C-SQUARED: ');
  108.    readln(m,n,c2);
  109.    IF ((n < m) OR (m < 0)) THEN BEGIN
  110.       writeln('Improper arguments');
  111.       GOTO 1
  112.    END;
  113.    factr := 1.0;
  114.    IF (m <> 0) THEN BEGIN
  115.       q1 := n;
  116.       FOR i := 1 to m DO BEGIN
  117.          factr := -0.5*factr*(n+i)*(q1/i);
  118.          q1 := q1-1.0
  119.       END
  120.    END;
  121.    v1[1] := n*(n+1)-m*(m+1)+c2/2.0;
  122.    IF (((n-m) MOD 2) = 0) THEN v2[1] := factr
  123.    ELSE v2[1] := -factr;
  124.    v2[2] := v1[1];
  125.    delv1[1] := delta*v1[1];
  126.    delv2[1] := delta*factr;
  127.    delv2[2] := delv1[1];
  128.    h1 := 0.1;
  129.    hmin := 0.0;
  130.    x1 := -1.0+dx;
  131.    x2 := 1.0-dx;
  132.    xf := 0.0;
  133.    writeln;
  134.    writeln('mu(-1)':26,'y(1-dx)':20,'mu(+1)':19);
  135.    REPEAT
  136.       shootf(nvar,v1,v2,delv1,delv2,n1,n2,x1,x2,
  137.          xf,eps,h1,hmin,f,dv1,dv2);
  138.       writeln;
  139.       writeln('v ':6,v1[1]:20:6,v2[1]:20:6,v2[2]:20:6);
  140.       writeln('dv':6,dv1[1]:20:6,dv2[1]:20:6,dv2[2]:20:6);
  141.    UNTIL (abs(dv1[1]) <= abs(eps*v1[1]))
  142. END.
  143.