home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / SFROID.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1991-04-29  |  2.3 KB  |  92 lines
  1. PROGRAM sfroid(input,output);
  2. LABEL 99;
  3. CONST
  4.    ne=3;
  5.    m=41;
  6.    nb=1;
  7.    nsi=ne;
  8.    nyj=ne;
  9.    nyk=m;
  10.    nci=ne;
  11.    ncj=3;      (* ncj=ne-nb+1 *)
  12.    nck=42;      (* nck=m+1 *)
  13.    nsj=7;      (* nsj=2*ne+1 *)
  14. TYPE
  15.    glyarray = ARRAY [1..ne,1..m] OF real;
  16.    glcarray = ARRAY [1..nci,1..ncj,1..nck] OF real;
  17.    glsarray = ARRAY [1..nsi,1..nsj] OF real;
  18.    glscalv = ARRAY [1..ne] OF real;
  19.    glindex = ARRAY [1..ne] OF integer;
  20. VAR
  21.    mm,n: integer;
  22.    i,itmax,k: integer;
  23.    h,c2,anorm: real;
  24.    conv,deriv,fac1,fac2: real;
  25.    q1,slowc: real;
  26.    scalv: glscalv;
  27.    indexv: glindex;
  28.    y: glyarray;
  29.    c: glcarray;
  30.    s: glsarray;
  31.    x: ARRAY [1..m] OF real;
  32. (*$I MODFILE.PAS *)
  33. (*$I PLGNDR.PAS *)
  34. (*$I DIFEQ.PAS *)
  35. (*$I RED.PAS *)
  36. (*$I PINVS.PAS *)
  37. (*$I BKSUB.PAS *)
  38. (*$I SOLVDE.PAS *)
  39. BEGIN
  40.    itmax := 100;
  41.    c2 := 0.0;
  42.    conv := 5.0e-6;
  43.    slowc := 1.0;
  44.    h := 1.0/(m-1);
  45.    writeln('Enter m n');
  46.    readln(mm,n);
  47.    IF (((n+mm) MOD 2) = 1) THEN BEGIN
  48.       indexv[1] := 1;
  49.       indexv[2] := 2;
  50.       indexv[3] := 3
  51.    END ELSE BEGIN
  52.       indexv[1] := 2;
  53.       indexv[2] := 1;
  54.       indexv[3] := 3
  55.    END;
  56.    anorm := 1.0;
  57.    IF (mm <> 0) THEN BEGIN
  58.       q1 := n;
  59.       FOR i := 1 TO mm DO BEGIN
  60.          anorm := -0.5*anorm*(n+i)*(q1/i);
  61.          q1 := q1-1.0
  62.       END
  63.    END;
  64.    FOR k := 1 TO (m-1) DO BEGIN
  65.       x[k] := (k-1)*h;
  66.       fac1 := 1.0-sqr(x[k]);
  67.       fac2 := exp((-mm/2.0)*ln(fac1));
  68.       y[1,k] := plgndr(n,mm,x[k])*fac2;
  69.       deriv := -((n-mm+1)*plgndr(n+1,mm,x[k])-
  70.          (n+1)*x[k]*plgndr(n,mm,x[k]))/fac1;
  71.       y[2,k] := mm*x[k]*y[1,k]/fac1+deriv*fac2;
  72.       y[3,k] := n*(n+1)-mm*(mm+1)
  73.    END;
  74.    x[m] := 1.0;
  75.    y[1,m] := anorm;
  76.    y[3,m] := n*(n+1)-mm*(mm+1);
  77.    y[2,m] := (y[3,m]-c2)*y[1,m]/(2.0*(mm+1.0));
  78.    scalv[1] := abs(anorm);
  79.    IF (y[2,m] > abs(anorm)) THEN scalv[2] := y[2,m] ELSE scalv[2] := abs(anorm);
  80.    IF (y[3,m] > 1.0) THEN scalv[3] := y[3,m] ELSE scalv[3] := 1.0;
  81.    WHILE true DO BEGIN
  82.       writeln('Enter c**2 or 999 to end.');
  83.       readln(c2);
  84.       IF (c2 = 999) THEN GOTO 99;
  85.       solvde(itmax,conv,slowc,scalv,indexv,ne,nb,m,y,nyj,nyk,
  86.          c,nci,ncj,nck,s,nsi,nsj);
  87.       writeln;
  88.       writeln('m = ',mm:2,'  n = ',n:2,'  c**2 = ',c2:7:3,
  89.          ' lam = ',y[3,1]+mm*(mm+1):10:6);
  90.    END;
  91. 99:   END.
  92.