c Based on algorithms from the book on c Numerical Recipes by Press et al. c All of these routines differ from the actual code used in c that book c and are unlikely to work as the c original routines did. c KENDALL'S NON-PARAMETRIC TEST OF CORRELATION c based on algorithm in Numerical Recipes P. 493 subroutine kendl(xfit, yfit, n) include 'typecom' include 'robcom' double precision xfit(*), yfit(*) integer n integer rn1, rn2, j, k c square root of two parameter (srt = 1.4142136d0) if(n .le. 1)then call xtext('NOT ENOUGH POINTS TO CALCULATE STATISTIC!') write(tstring, *)'Number of points =', n call xtext(tstring) return end if rn1 = 0 rn2 = 0 is = 0 do 120 j = 1, n - 1 do 110 k = j + 1, n aa1 = xfit(j) - xfit(k) aa2 = yfit(j) - yfit(k) a12 = aa1 * aa2 if(a12 .ne. 0.0d0)then rn1 = rn1 + 1 rn2 = rn2 + 1 if(a12 .gt. 0.0d0)then is = is + 1 else is = is - 1 end if else if(aa1 .ne. 0.0d0) rn1 = rn1 + 1 if(aa2 .ne. 0.0d0) rn2 = rn2 + 1 end if 110 continue 120 continue t = dble(is) / sqrt(dble(rn1) * dble(rn2)) var = (4.0d0 * n + 10.0d0) / (9.0d0 * n * (n - 1.0d0)) z = t / sqrt(var) pb = erfrc(abs(z)/srt) write(tstring,*)'Kendall tau statistic = ',t call totext(tstring) call write8(tstring) write(tstring,*)'Probability = ',pb call totext(tstring) call write8(tstring) end C ERROR FUNCTION c based on algorithm in Numerical Recipes P.164 double precision function erfrc(xin) include 'typecom' double precision xin parameter (HALF = 0.5d0) y = abs(xin) t1 = 1.0d0 / (1.0d0 + HALF * y) erfrc = t1*exp(-y*y - 1.26551223d0 + t1 * + (1.00002368d0 + t1 * (0.37409196d0 + + t1 * (0.09678418d0 + t1 * (-0.18628806d0 + t1 * + (0.27886807d0 + t1 * (-1.13520398d0 + + t1 * (1.48851587d0 + t1 * (-0.82215223d0 + + t1 *0.17087277d0))))))))) if(xin .lt. 0.0d0) erfrc = 2.0d0 - erfrc end C MEDIAN OF DATA ARRAY C based on algorithm in Numerical Recipes P. 461 c added kludge of escape route to avoid infinite loops subroutine mdian(x, guess, n) include 'typecom' include 'robcom' double precision x(*) c guess is supplied by the calling routine as first guess for median double precision guess integer n parameter(HUGE = 1.0d30, afct = 1.5d0, amp = 1.5d0) parameter(MAXLOOPS = 1000) parameter (HALF = 0.5d0) integer itimes double precision sumx, sum, dum, xx, a, esp c a = HALF * (x(1) + x(n)) a = guess esp = abs(x(n) - x(1)) ap = HUGE am = -HUGE itimes = 0 10 continue c print*,'itimes = ', itimes itimes = itimes + 1 c print*,'0.5 * (xp + xm) = ', 0.5d0 * (xp + xm) sum = 0.0d0 sumx = 0.0d0 np = 0 nm = 0 xp = HUGE xm = -HUGE do 110 j = 1, n xx = x(j) if(xx .ne. a)then if(xx .gt. a)then np = np + 1 if(xx .lt. xp)xp = xx else if(xx .lt. a)then nm = nm + 1 if(xx .gt. xm)xm = xx end if dum = 1.0d0 / (esp + abs(xx - a)) sum = sum + dum sumx = sumx + xx * dum end if 110 continue c guess too low if((np - nm) .ge. 2 .and. itimes .lt. MAXLOOPS)then am = a aa = xp + max(0.0d0, sumx/sum-a) * amp if(aa .gt. ap)aa = HALF*(a + ap) esp = afct * abs(aa - a) a = aa goto 10 c guess too high else if((nm - np) .ge. 2 .and. itimes .lt. MAXLOOPS)then ap = a aa = xm + min(0.0d0, sumx/sum-a) * amp if(aa .lt. am)aa = HALF * (a + am) esp = afct * abs(aa - a) a = aa goto 10 else if(mod(n, 2) .eq. 0)then if(np.eq.nm)then xmed = HALF * (xp + xm) else if(np .gt. nm)then xmed = HALF * (a + xp) else xmed = HALF*(xm + a) end if else if(np .eq. nm)then xmed = a else if(np .gt. nm)then xmed = xp else xmed = xm end if end if end if if(itimes .ge. MAXLOOPS)then call totext('WARNING: Exiting infinite loop') call totext('MEDIAN MAY HAVE INCORRECT VALUE') end if write(tstring,*)'Median = ',xmed call totext(tstring) call write8(tstring) end C GAUSSIAN RANDOM NUMBER GENERATOR C based on algorithm in Numerical Recipes P. 203 double precision function gasd(dummy) include 'typecom' include 'robcom' integer dummy data ist/0/ if(ist .eq. 0)then 10 val1 = 2.0d0 * ran1(dummy) - 1.0d0 val2 = 2.0d0 * ran1(dummy) - 1.0d0 rad = val1 * val1 + val2 * val2 if(rad .ge. 1.0d0)goto 10 fact = sqrt(-2.0d0 * log(rad) / rad) gst = val1 * fact gasd = val2 * fact ist = 1 else gasd = gst ist = 0 end if end c Numerical recipes algorithmq didn't work (dunno why) double precision function ran1(dummy) include 'typecom' integer*4 dummy ran1 = rann(dummy) end C RAN1: RETURNS UNIFORM RANDOM DEVIATE BETWEEN 0 AND 1 C based on algorithm in Numerical Recipes P.196 double precision function rann(dummy) include 'typecom' integer dummy include 'robcom' double precision xa(97) integer m1, m2, m3 integer ia1, ia2, ia3 double precision rm1, rm2 parameter(m1 = 259200, ia1 = 7141, ik1 = 54773, rm1 = 1.0d0/m1) parameter(m2 = 134456, ia2 = 8121, ik2 = 28411, rm2 = 1.0d0/m2) parameter(m3 = 243000, ia3 = 4561, ik3 = 51349) data ifrc/0/ if(dummy .lt. 0 .or. ifrc .eq. 0)then ifrc = 1 jx1 = mod(ik1-dummy,m1) jx1 = mod(ia1*jx1+ik1,m1) jx2 = mod(jx1,m2) jx1 = mod(ia1*jx1+ik1,m1) jx3 = mod(jx1,m3) do 110 j = 1,97 jx1 = mod(ia1*jx1+ik1,m1) jx2 = mod(ia2*jx2+ik2,m2) xa(j) = (dble(jx1) + dble(jx2) * rm2) * rm1 110 continue dummy = 1 end if jx1 = mod(ia1 * jx1 + ik1,m1) jx2 = mod(ia2 * jx2 + ik2,m2) jx3 = mod(ia3 * jx3 + ik3,m3) j = 1 + (97 * jx3) / m3 rann = xa(j) xa(j) = (dble(jx1) + dble(jx2) * rm2) * rm1 end C CUBIC SPLINE ROUTINE c based on algorithm in Numerical Recipes P.88 subroutine nrspln(x, y, n, yp1, ypn, yout, u) include 'typecom' include 'robcom' double precision x(*), y(*), yout(*), u(*) parameter (HUGE = 0.99d+30) parameter (HALF = 0.5d0) if(yp1 .gt. HUGE)then yout(1) = 0.0d0 u(1) = 0.0d0 else yout(1) = -HALF u(1) = (3.0d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) end if do 110 i = 2, n - 1 sig = (x(i)-x(i-1))/(x(i+1)-x(i-1)) p = sig*yout(i-1) + 2.0d0 yout(i) = (sig-1.0d0)/p u(i) = (6.0d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) +/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p 110 continue if(ypn .gt. HUGE)then pn = 0.0d0 un = 0.0d0 else pn = HALF un = (3.0d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) end if yout(n) = (un - pn * u(n - 1))/(pn * yout(n - 1) + 1.0d0) do 120 index = n - 1, 1, -1 yout(index) = yout(index) * yout(index + 1) + u(index) 120 continue end C INTERPOLATE RESULTS OF SPLINE FIT c based on algorithm in Numerical Recipes P.89 subroutine splnt(xin, yin, yin2, n, x, y) include 'typecom' include 'robcom' double precision xin(*), yin(*), yin2(*) integer lo, hi lo = 1 hi = n 10 if((hi - lo) .gt. 1)then kt = (hi + lo) / 2 if(xin(kt) .gt. x)then hi = kt else lo = kt end if goto 10 end if h = xin(hi) - xin(lo) if(h .eq. 0.0d0)then call xtext('ERROR: BAD INPUT DETECTED IN S/R SPLINT') return end if a = (xin(hi) - x) / h b = (x - xin(lo)) / h y = a * yin(lo) + b * yin(hi) + + ((a * a * a - a) * yin2(lo) + + (b * b * b - b) * yin2(hi)) * + (h * h)/6.0d0 end