home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / MRQMIN.DEM < prev    next >
Encoding:
Text File  |  1991-04-29  |  2.8 KB  |  116 lines
  1. PROGRAM D14R8(input,output);
  2. (* driver for routine MRQMIN *)
  3. LABEL 1;
  4. CONST
  5.    npt=100;
  6.    mma=6;
  7.    spread=0.001;
  8. TYPE
  9.    glmma = ARRAY [1..mma] OF real;
  10.    glnparam = glmma;
  11.    gllista = ARRAY [1..mma] OF integer;
  12.    glcovar = ARRAY [1..mma,1..mma] OF real;
  13.    glnalbynal = glcovar;
  14.    glncabynca = glcovar;
  15.    glnpbynp = glcovar;
  16.    glnpbymp = glcovar;
  17.    glndata = ARRAY [1..npt] OF real;
  18. VAR
  19.    glochisq : real;
  20.    glatry,glbeta : glmma;
  21.    glinext,glinextp : integer;
  22.    glma : ARRAY [1..55] OF real;
  23.    gliset : integer;
  24.    glgset : real;
  25.    alamda,chisq,ochisq : real;
  26.    i,idum,itst,j,jj,k,mfit : integer;
  27.    x,y,sig : glndata;
  28.    lista : gllista;
  29.    a,gues : glmma;
  30.    covar,alpha : glcovar;
  31.  
  32. (*$I MODFILE.PAS *)
  33. (*$I COVSRT.PAS *)
  34.  
  35. (*$I RAN3.PAS *)
  36.  
  37. (*$I GASDEV.PAS *)
  38.  
  39. PROCEDURE funcs(x: real; a: glnparam; VAR y: real;
  40.        VAR dyda: glnparam; na: integer);
  41. (* Programs using routine FGAUSS must define the type
  42. TYPE
  43.    glnparam = ARRAY [1..na] OF real;
  44. where na is the number of parameters. *)
  45. VAR
  46.    i,ii : integer;
  47.    fac,ex,arg : real;
  48. BEGIN
  49.    y := 0.0;
  50.    FOR ii := 1 to (na DIV 3) DO BEGIN
  51.       i := 3*ii-2;
  52.       arg := (x-a[i+1])/a[i+2];
  53.       ex := exp(-sqr(arg));
  54.       fac := a[i]*ex*2.0*arg;
  55.       y := y+a[i]*ex;
  56.       dyda[i] := ex;
  57.       dyda[i+1] := fac/a[i+2];
  58.       dyda[i+2] := fac*arg/a[i+2]
  59.    END
  60. END;
  61.  
  62. (*$I GAUSSJ.PAS *)
  63.  
  64. (*$I MRQCOF.PAS *)
  65.  
  66. (*$I MRQMIN.PAS *)
  67.  
  68. BEGIN
  69.    gliset := 0;
  70.    a[1] := 5.0; a[2] := 2.0; a[3] := 3.0;
  71.    a[4] := 2.0; a[5] := 5.0; a[6] := 3.0;
  72.    gues[1] := 4.5; gues[2] := 2.2; gues[3] := 2.8;
  73.    gues[4] := 2.5; gues[5] := 4.9; gues[6] := 2.8;
  74.    idum := -911;
  75.    FOR i := 1 to 100 DO BEGIN
  76.       x[i] := 0.1*i;
  77.       y[i] := 0.0;
  78.       FOR jj := 1 to 2 DO BEGIN
  79.          j := 3*jj-2;
  80.          y[i] := y[i]+a[j]*exp(-sqr((x[i]-a[j+1])/a[j+2]))
  81.       END;
  82.       y[i] := y[i]*(1.0+spread*gasdev(idum));
  83.       sig[i] := spread*y[i]
  84.    END;
  85.    mfit := 6;
  86.    FOR i := 1 to mfit DO lista[i] := i;
  87.    alamda := -1;
  88.    FOR i := 1 to mma DO a[i] := gues[i];
  89.    mrqmin(x,y,sig,npt,a,mma,lista,mfit,covar,alpha,
  90.       mma,chisq,alamda);
  91.    k := 1;
  92.    itst := 0;
  93. 1:   writeln;
  94.    writeln('Iteration #',k:2,'chi-squared:':17,chisq:10:4,
  95.       'alamda:':10,alamda:9);
  96.    writeln('a[1]':7,'a[2]':8,'a[3]':8,'a[4]':8,'a[5]':8,'a[6]':8);
  97.    FOR i := 1 to 6 DO write(a[i]:8:4);
  98.    writeln;
  99.    k := k+1;
  100.    ochisq := chisq;
  101.    mrqmin(x,y,sig,npt,a,mma,lista,mfit,covar,alpha,
  102.       mma,chisq,alamda);
  103.    IF (chisq > ochisq) THEN BEGIN
  104.       itst := 0
  105.    END ELSE IF (abs(ochisq-chisq) < 0.1) THEN BEGIN
  106.       itst := itst+1
  107.    END;
  108.    IF (itst < 2) THEN GOTO 1;
  109.    alamda := 0.0;
  110.    mrqmin(x,y,sig,npt,a,mma,lista,mfit,covar,alpha,
  111.       mma,chisq,alamda);
  112.    writeln('Uncertainties:');
  113.    FOR i := 1 to 6 DO write(sqrt(covar[i,i]):8:4);
  114.    writeln
  115. END.
  116.