- #1
NathanM
- 14
- 0
Hey guys! Having issues with a code I'm writing. It doesn't fit the typical format and I may have already accidentally posted it in the wrong forum, so here's a copy of what was there with the newest version of the code. It'll (eventually) be a four-parameter least squares minimisation, but for the moment I'm just trying to get the single parameter version to work. I've got the first module done but, when I run it it gets caught up in a loop and doesn't actually minimise anything! Good start huh?
Here's my code, bless me with your insight!
PROGRAM assignment
! A program designed to fit experiemental data, using the method
! of least squares to minimise the associated chi-squared and
! obtain the four control parameters A,B,h1 and h2.
!*****************************************************************
IMPLICIT NONE
INTEGER i,t0
DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
DOUBLE PRECISION deltah,Cs(17)
OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
DO i=1,17
READ(21,*)t(i),Ct(i),eCt(i)
END DO
CLOSE(21)
!Read in data.txt as three one dimensional arrays.
!*****************************************************************
!OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
!DO i=1,17
! WRITE(21,*)t(i),Ct(i),eCt(i)
!END DO
!CLOSE(21)
!
!Just to check input file is being read correctly.
!*****************************************************************
!****************************Parameters***************************
Aa= 0
Ba= 0
h1a= 0
h2a= 0
!*****************************************************************
!**********************Minimising Lamda1 (h1)*********************
deltah= 0.0001
h1b= 0.001
h1a= 0
DO 10
chisqa= 0
DO 20 i= 1,17
Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
Ctdiff= Ct - Cs
c= Ctdiff**2/eCt**2
chisqa= chisqa + c(i)
20 END DO
!Use initial h1 value of 0 to calculate start-point chisq.
chisqb= 0
DO 30 i= 1,17
h1b= h1b + deltah
Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
Ctdiff(i)= Ct(i) - Cs(i)
c= Ctdiff(i)**2/eCt(i)**2
chisqb= chisqb + c(i)
30 END DO
!First-step h1 used to find competing chisq for comparison.
WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1a
WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1b
!Prints the two calculated chisq values to screen.
dchisq= chisqa - chisqb
IF (dchisq.GT.0) THEN
h1a= h1b
ELSE IF (dchisq.LE.0) THEN
deltah= deltah-((deltah*2)/100)
END IF
IF (chisqb.LE.6618.681) EXIT
10 END DO
WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
!*****************************************************************
END PROGRAM assignmentThanks for the help.
Nathan
EDIT: I should be getting a chi-squared of 6618.681 from this, but it's just stuck between 6921.866 and 6920.031. Help!
Here's my code, bless me with your insight!
PROGRAM assignment
! A program designed to fit experiemental data, using the method
! of least squares to minimise the associated chi-squared and
! obtain the four control parameters A,B,h1 and h2.
!*****************************************************************
IMPLICIT NONE
INTEGER i,t0
DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
DOUBLE PRECISION deltah,Cs(17)
OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
DO i=1,17
READ(21,*)t(i),Ct(i),eCt(i)
END DO
CLOSE(21)
!Read in data.txt as three one dimensional arrays.
!*****************************************************************
!OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
!DO i=1,17
! WRITE(21,*)t(i),Ct(i),eCt(i)
!END DO
!CLOSE(21)
!
!Just to check input file is being read correctly.
!*****************************************************************
!****************************Parameters***************************
Aa= 0
Ba= 0
h1a= 0
h2a= 0
!*****************************************************************
!**********************Minimising Lamda1 (h1)*********************
deltah= 0.0001
h1b= 0.001
h1a= 0
DO 10
chisqa= 0
DO 20 i= 1,17
Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
Ctdiff= Ct - Cs
c= Ctdiff**2/eCt**2
chisqa= chisqa + c(i)
20 END DO
!Use initial h1 value of 0 to calculate start-point chisq.
chisqb= 0
DO 30 i= 1,17
h1b= h1b + deltah
Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
Ctdiff(i)= Ct(i) - Cs(i)
c= Ctdiff(i)**2/eCt(i)**2
chisqb= chisqb + c(i)
30 END DO
!First-step h1 used to find competing chisq for comparison.
WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1a
WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1b
!Prints the two calculated chisq values to screen.
dchisq= chisqa - chisqb
IF (dchisq.GT.0) THEN
h1a= h1b
ELSE IF (dchisq.LE.0) THEN
deltah= deltah-((deltah*2)/100)
END IF
IF (chisqb.LE.6618.681) EXIT
10 END DO
WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
!*****************************************************************
END PROGRAM assignmentThanks for the help.
Nathan
EDIT: I should be getting a chi-squared of 6618.681 from this, but it's just stuck between 6921.866 and 6920.031. Help!