- #1

sarah_frtrn

- 6

- 1

Fortran:

```
program main
!change V, B, and P in ROOT_CHISQ subroutine to see what answers
implicit none
real :: A, B, V, final
integer :: N
a=0
b=12.5
v=10
n=1000000
call simp_p(A,B,N,V,final)
print *, final
call root_chisq()
end program main
subroutine simp_p(A,B,N,V,final)
real :: A,B ! endpoints and number of divisions
real :: DX, integral, x, place
real :: V ! degrees freedom
real :: answer, final ! return from subroutine
integer :: N, i
integer(selected_int_kind(r=12)) :: i12
! calculate partition width
DX = (B-A)/FLOAT(N)
!print *, "DX = ", DX
!print *, "V = ", V
! initialize integral variable
integral = 0
call equation(A, V, answer)
!print *, "answer a = ", answer
if (answer.ne.answer) THEN
integral = integral + 0
else if (answer.ge.huge(i12)) THEN
integral = integral + 0
else
integral = integral + answer
end if
call equation(B, V, answer)
!print *, "answer b = ", answer
if (answer.ne.answer) THEN
integral = integral + 0
else if (answer.ge.huge(i12)) THEN
integral = integral + 0
else
integral = integral + answer
end if
!now compute values of remaining partitions
do i=1, N-1
!print *, "i equals = ", i
x = A+i*DX
!print *, "xval = ", x
call equation(x, v, answer )
if (answer.ne.answer) THEN
cycle
else if (answer.ge.huge(i12)) THEN
cycle
end if
if (MOD(i,2).eq.0) THEN
integral = integral + 2*answer
else if (mod(i,2).ne.0) THEN
INTEGRAL = INTEGRAL + 4*answer
end if
!print *, integral
!print *,
!print *,
end do
!print *, 'final: ', integral
integral = integral * DX/3
! EVERYTHING BELOW HAS BEEN CHECKED, SOMETHING IS WRONG
! WITH THE INTEGRATION
call calc_gamma(v/2, answer)
place = (2**(v/2))*answer
place = place**(-1)
!print *, "place = ", place
final= place*integral
print *, 'final ANSWER = ', final
end subroutine simp_p
subroutine equation(x, v, answer)
real, intent(in) :: v, x
real, intent(out) :: answer
real :: pt1, pt2
pt1 = x**((v/2)-1)
pt2 = exp(-x/2)
answer = pt1*pt2
!print *, "answer from eq = ", answer
return
end subroutine equation
subroutine calc_gamma(x, answer)
real :: x, answer, difference, previous_answer
integer :: n, i
answer = 1
n = 1
difference = 1
previous_answer = 1
do i=1, 56140
answer = answer * (((1/float(n)) + 1)**x)/((x/float(n))+ 1)
difference = abs(previous_answer - answer)
n=n+1
end do
answer = answer * (1/x)
end subroutine calc_gamma
subroutine root_chisq()
! assuming root lies somewhere in the interval 0,26
real :: p = 0.75
real :: v=10
real :: error = 0.01
real :: a=0, answera, answerb, m, answerm
real :: b=100.0 ! highest possible cutoff value when v=10
real :: answer
integer :: N=1000
logical :: check=.false.
!first check value of function at two endpoints
call simp_p(a,a,N,V, answer)
print *, "answera = ", answer
print *, "answer - p = ", answer-p
answera = answer-p
if (abs(answera).LE.0.0001) THEN
print *, "ANSWER FOUND"
else
print *, "KEEP SEARCHING"
end if
call simp_p(a,b,n,v,answer)
answerb = answer-p
print *, "answer = ", answer
if (abs(answerb).LE.0.0001) THEN
print *, "ANSWER FOUND"
check=.true.
else
print *, "KEEP SEARCHING"
end if
do
m=(a+b)/2
call simp_p(a,m,n,v,answer)
answerm = answer-p
call simp_p(a,a,n,v,answer)
answera = answer-p
call simp_p(a,b,n,v,answer)
answerb = answer-p
print *, "answer a = ", answera
print *, "answer b = ", answerb
print *, "answer m = ", answerm
if(abs(answerm).le.0.001) THEN
print *, "FINAL ANSWERM = ", answerm
exit
else if (abs(answerb).le.0.001) THEN
print *, "FINAL ANSWERB = ", answerb
exit
else if (abs(answera).le.0.001) THEN
print *, "FINAL ANSWERA = ", answera
exit
else
if ((answera.lt.0).AND.(answerm.GT.0)) THEN
b=m
answerb=answerm
print *, "replacing b"
print *, "m = ", (a+b)/2
else
a=m
answera=answerm
print *, "replacing a"
print *, "m = ", (a+b)/2
end if
end if
if (abs(b-a).LT.e) THEN
print *, "final answerb = ", b
print *, "final answera = ", a
print *, "answerb = ", answerb
print *, "answera = ", answera
print *, (a+b)/2
print *, b
exit
end if
end do
return
end subroutine root_chisq
```