program gauss_tst integer array parameter(array=200) integer i,j,nmbr,rot,rotc real b1(array),b1ave,b1adev,b1sdev,b1var,b1skew,b1curt real b2(array),b2ave,b2adev,b2sdev,b2var,b2skew,b2curt real rho(array),rave,radev,rsdev,rvar,rskew,rcurt real a1(array),a2(array),sig(array),l(array) real cor12,cov12,cor13,cov13,cor23,cov23 real prob12,test12,prob13,test13,prob23,test23 real cor(3,3),e(3),ev(3,3) real cov(3,3),ec(3),evc(3,3) real tan integer ha1(9),ha2(9),ha1_left,ha2_left read(5,*) nmbr do i=1,nmbr read(5,*) tan,a1(i),a2(i),sig(i),l(i) enddo call data_form(nmbr,a1,a2,l,b1,b2,rho) call moment(b1,nmbr,b1ave,b1adev,b1sdev,b1var,b1skew,b1curt) call moment(b2,nmbr,b2ave,b2adev,b2sdev,b2var,b2skew,b2curt) call moment(rho,nmbr,rave,radev,rsdev,rvar,rskew,rcurt) c call histogram(nmbr,b1,b2,b1ave,b2ave,b1adev,b2adev c & ,ha1,ha2,ha1_left,ha2_left) call pearsn(b1,b2,nmbr, cor12, prob12, test12) call pearsn(b1,rho,nmbr, cor13, prob13, test13) call pearsn(b2,rho,nmbr, cor23, prob23, test23) cov12=cor12*b1sdev*b2sdev cov13=cor13*b1sdev*rsdev cov23=cor23*b2sdev*rsdev cov(1,1)=b1var cov(1,2)=cov12 cov(1,3)=cov13 cov(2,1)=cov(1,2) cov(2,2)=b2var cov(2,3)=cov23 cov(3,1)=cov(1,3) cov(3,2)=cov(2,3) cov(3,3)=rvar cor(1,1)=1.0 cor(1,2)=cor12 cor(1,3)=cor13 cor(2,1)=cor(1,2) cor(2,2)=1.0 cor(2,3)=cor23 cor(3,1)=cor(1,3) cor(3,2)=cor(2,3) cor(3,3)=1.0 call jacobi(cor,3,3,e,ev,rot) call jacobi(cov,3,3,ec,evc,rotc) write(6,*) 'i,e(i),ev(1,i),ev(2,i),ev(3,i)' write(6,*) '1',e(1),ev(1,1),ev(2,1),ev(3,1) write(6,*) '2',e(2),ev(1,2),ev(2,2),ev(3,2) write(6,*) '3',e(3),ev(1,3),ev(2,3),ev(3,3) call eigsrt(e,ev,3,3) write(6,*) 'now sorted' write(6,*) 'i,e(i),ev(1,i),ev(2,i),ev(3,i)' write(6,*) '1',e(1),ev(1,1),ev(2,1),ev(3,1) write(6,*) '2',e(2),ev(1,2),ev(2,2),ev(3,2) write(6,*) '3',e(3),ev(1,3),ev(2,3),ev(3,3) write(6,*) 'i,ec(i),evc(1,i),evc(2,i),evc(3,i)' write(6,*) '1',ec(1),evc(1,1),evc(2,1),evc(3,1) write(6,*) '2',ec(2),evc(1,2),evc(2,2),evc(3,2) write(6,*) '3',ec(3),evc(1,3),evc(2,3),evc(3,3) write(6,*) 'cor12,cov12,prob12 = ',cor12,cov12,prob12 write(6,*) 'cor13,cov13,prob13 = ',cor13,cov13,prob13 write(6,*) 'cor23,cov23,prob23 = ',cor23,cov23,prob23 write(6,*) 'i,ha1(i),ha2(i)' do i=1,9 write(6,*) i,ha1(i),ha2(i) enddo write(6,*) 'ha1_left,ha2_left = ',ha1_left,ha2_left write(6,*) 'nmbr = ',nmbr write(6,*) 'b1ave = ',b1ave write(6,*) 'b1adev = ',b1adev write(6,*) 'b1sdev = ',b1sdev write(6,*) 'b1var = ',b1var write(6,*) 'b1skew = ',b1skew write(6,*) 'b1curt = ',b1curt write(6,*) 'b2ave = ',b2ave write(6,*) 'b2adev = ',b2adev write(6,*) 'b2sdev = ',b2sdev write(6,*) 'b2var = ',b2var write(6,*) 'b2skew = ',b2skew write(6,*) 'b2curt = ',b2curt write(6,*) 'rave = ',rave write(6,*) 'radev = ',radev write(6,*) 'rsdev = ',rsdev write(6,*) 'rvar = ',rvar write(6,*) 'rskew = ',rskew write(6,*) 'rcurt = ',rcurt do i=1,nmbr write(6,*) b1(i),b2(i),rho(i) enddo stop end c************************************************************************ subroutine data_form(n,c1,c2,ell,d1,d2,p) c************************************************************************ integer array parameter(array=200) integer n,i real c1(array),c2(array),d1(array),d2(array) real ell(array),p(array) real zbrent,roott external roott real a1,a2,ell_0 common /state/ a1,a2,ell_0 do i=1,n d1(i)=c1(i) d2(i)=sqrt(abs(c2(i))) ell_0=ell(i) a1=c1(i) a2=c2(i) p(i)=zbrent(roott,1.0,ell_0,0.01) c p(i)=sqrt(p(i)) enddo return end c************************************************************************ real function roott(x) c************************************************************************ real x real a1,a2,ell_0 real bound complex dl,dz common /state/ a1,a2,ell_0 common /int_data/ bound dz=a1+2.0*a2*x dl=sqrt(1.0+dz**2) if (a2 .eq. 0.0) then roott=ell_0-dl*x else roott=ell_0-((dz*sqrt(1+dz**2))-a1*sqrt(1+a1**2) & -asinh(a1)+asinh(dz)) & /(4.0*a2) endif C bound=0.02 C call integral(dl,0.0,x,ellc) C C ell=real(ellc) C C roott=ell_0-ell return end Cc************************************************************************ C C complex function dl(xx) C Cc************************************************************************ C C real xx C real a1,a2,ell_0 C real dz,dll C C C common /state/ a1,a2,ell_0 C C dz=a1+0.9*a2*xx**(-0.1) C C dll=sqrt(1.0+dz**2) C C dl=cmplx(dll,0.0) C C return C end c************************************************************************ real function asinh(x) c************************************************************************ real x asinh=alog(x + sqrt(x**2+1.0)) return end c************************************************************************ subroutine histogram(nmbr,b1,b2,b1ave,b2ave,b1adev,b2adev & ,ha1,ha2,ha1_left,ha2_left) c************************************************************************ integer array parameter(array=200) integer j,nmbr integer ha1(9),ha2(9),ha1_left,ha2_left real b1(array),b2(array),b1ave,b2ave,b1adev,b2adev do j=1,9 ha1(j)=0 ha2(j)=0 ha1_left=0 ha2_left=0 enddo do j=1,nmbr if((b1(j) .lt. b1ave + b1adev/4).and. & (b1(j) .gt. b1ave - b1adev/4))then ha1(5)=ha1(5)+1 else if((b1(j) .lt. b1ave - b1adev/4.0).and. & (b1(j) .gt. b1ave - 3.0*b1adev/4.0))then ha1(4)=ha1(4)+1 else if((b1(j) .lt. b1ave - 3.0*b1adev/4.0).and. & (b1(j) .gt. b1ave - 5.0*b1adev/4.0))then ha1(3)=ha1(3)+1 else if((b1(j) .lt. b1ave - 5.0*b1adev/4.0).and. & (b1(j) .gt. b1ave - 7.0*b1adev/4.0))then ha1(2)=ha1(2)+1 else if((b1(j) .lt. b1ave - 7.0*b1adev/4.0).and. & (b1(j) .gt. b1ave - 9.0*b1adev/4.0))then ha1(1)=ha1(1)+1 else if((b1(j) .gt. b1ave + b1adev/4.0).and. & (b1(j) .lt. b1ave + 3.0*b1adev/4.0))then ha1(6)=ha1(6)+1 else if((b1(j) .gt. b1ave + 3.0*b1adev/4.0).and. & (b1(j) .lt. b1ave + 5.0*b1adev/4.0))then ha1(7)=ha1(7)+1 else if((b1(j) .gt. b1ave + 5.0*b1adev/4.0).and. & (b1(j) .lt. b1ave + 7.0*b1adev/4.0))then ha1(8)=ha1(8)+1 else if((b1(j) .gt. b1ave + 7.0*b1adev/4.0).and. & (b1(j) .lt. b1ave + 9.0*b1adev/4.0))then ha1(9)=ha1(9)+1 else ha1_left=ha1_left+1 endif if((b2(j) .lt. b2ave + b2adev/4).and. & (b2(j) .gt. b2ave - b2adev/4))then ha2(5)=ha2(5)+1 else if((b2(j) .lt. b2ave - b2adev/4.0).and. & (b2(j) .gt. b2ave - 3.0*b2adev/4.0))then ha2(4)=ha2(4)+1 else if((b2(j) .lt. b2ave - 3.0*b2adev/4.0).and. & (b2(j) .gt. b2ave - 5.0*b2adev/4.0))then ha2(3)=ha2(3)+1 else if((b2(j) .lt. b2ave - 5.0*b2adev/4.0).and. & (b2(j) .gt. b2ave - 7.0*b2adev/4.0))then ha2(2)=ha2(2)+1 else if((b2(j) .lt. b2ave - 7.0*b2adev/4.0).and. & (b2(j) .gt. b2ave - 9.0*b2adev/4.0))then ha2(1)=ha2(1)+1 else if((b2(j) .gt. b2ave + b2adev/4.0).and. & (b2(j) .lt. b2ave + 3.0*b2adev/4.0))then ha2(6)=ha2(6)+1 else if((b2(j) .gt. b2ave + 3.0*b2adev/4.0).and. & (b2(j) .lt. b2ave + 5.0*b2adev/4.0))then ha2(7)=ha2(7)+1 else if((b2(j) .gt. b2ave + 5.0*b2adev/4.0).and. & (b2(j) .lt. b2ave + 7.0*b2adev/4.0))then ha2(8)=ha2(8)+1 else if((b2(j) .gt. b2ave + 7.0*b2adev/4.0).and. & (b2(j) .lt. b2ave + 9.0*b2adev/4.0))then ha2(9)=ha2(9)+1 else ha2_left=ha2_left+1 endif enddo return end SUBROUTINE EIGSRT(D,V,N,NP) DIMENSION D(NP),V(NP,NP) DO 13 I=1,N-1 K=I P=D(I) DO 11 J=I+1,N IF(D(J).GE.P)THEN K=J P=D(J) ENDIF 11 CONTINUE IF(K.NE.I)THEN D(K)=D(I) D(I)=P DO 12 J=1,N P=V(J,I) V(J,I)=V(J,K) V(J,K)=P 12 CONTINUE ENDIF 13 CONTINUE RETURN END