home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / SPEAR.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1991-04-29  |  1000 b   |  35 lines
  1. PROCEDURE spear(data1,data2: narray; n: integer;
  2.       VAR wksp1,wksp2: narray;
  3.       VAR d,zd,probd,rs,probrs: real);
  4. (* Programs using routine SPEAR must define types
  5. TYPE
  6.    narray = ARRAY [1..n] OF real;
  7.    glsarray = narray;
  8. in the calling routine *)
  9. VAR
  10.    j: integer;
  11.    vard,t,sg,sf,fac,en3n,en,df,aved: real;
  12. BEGIN
  13.    FOR j := 1 TO n DO BEGIN
  14.       wksp1[j] := data1[j];
  15.       wksp2[j] := data2[j]
  16.    END;
  17.    sort2(n,wksp1,wksp2);
  18.    crank(n,wksp1,sf);
  19.    sort2(n,wksp2,wksp1);
  20.    crank(n,wksp2,sg);
  21.    d := 0.0;
  22.    FOR j := 1 TO n DO d := d+sqr(wksp1[j]-wksp2[j]);
  23.    en := n;
  24.    en3n := en*en*en-en;
  25.    aved := en3n/6.0-(sf+sg)/12.0;
  26.    fac := (1.0-sf/en3n)*(1.0-sg/en3n);
  27.    vard := ((en-1.0)*sqr(en)*sqr(en+1.0)/36.0)*fac;
  28.    zd := (d-aved)/sqrt(vard);
  29.    probd := erfcc(abs(zd)/1.4142136);
  30.    rs := (1.0-(6.0/en3n)*(d+0.5*(sf+sg)))/fac;
  31.    t := rs*sqrt((en-2.0)/((1.0+rs)*(1.0-rs)));
  32.    df := en-2.0;
  33.    probrs := betai(0.5*df,0.5,df/(df+sqr(t)))
  34. END;
  35.