PROGRAM FTCTempAnalysis3

! Written by Kate T-C
! August 27, 2009
! This program calculates a trend and correlation coefficient in 
! a 1-dimensional array of data

implicit none

integer, parameter :: years = 41, months=12
real, dimension(years*months) :: TimeIndexes, AvgMonthlyTemps
real, dimension(years) :: YearIndexes, JanTemps, JulTemps
integer :: counter = 0, n=years*months
real :: slope = 0, intercept = 0, rsquare=0
real :: CalcCorrCoef

    ! Read in my data from the command line
	do counter=1,years*months
		read *, AvgMonthlyTemps(counter)  
		TimeIndexes(counter) = counter
	end do

	! Pull out January's and July's
	JanTemps = AvgMonthlyTemps(1:years*months:12)
	JulTemps = AvgMonthlyTemps(7:years*months:12)
	YearIndexes = floor(TimeIndexes(1:years*months:12)/12.0)+1.0
	
	! Print statistics for January
	call CalcSlopeIntercept(years, YearIndexes, JanTemps, slope, intercept)
	rsquare = CalcCorrCoef(years, YearIndexes, JanTemps, slope, intercept)
	print *, 'For January:'
	print *, 'slope = ', slope, ' intercept is ', intercept, ' with r**2= ', rsquare
	
	! Print statistics for July
	call CalcSlopeIntercept(years, YearIndexes, JulTemps, slope, intercept)
	rsquare = CalcCorrCoef(years, YearIndexes, JulTemps, slope, intercept)
	print *, 'For July:'
	print *, 'slope = ', slope, ' intercept is ', intercept, ' with r**2= ', rsquare
	
END

SUBROUTINE CalcSlopeIntercept(n,xSeries,yValues,slope,intercept)

implicit none

integer, intent(IN) :: n
real, intent(IN), dimension(n) :: xSeries, yValues
real, intent(OUT) :: slope, intercept
integer :: counter = 0
real :: SumX=0, SumXSquare=0, SumY=0, SumXY=0

	! Calculate the four sums that we need
	SumX = sum(xSeries)
	SumY = sum(yValues)
		
	SumXSquare = 0
	SumXY = 0
	do counter=1,n
		SumXSquare = SumXSquare + xSeries(counter)**2
		SumXY = SumXY + yValues(counter)*xSeries(counter)
	end do
	
	! Calculate the slope of the trend
	slope = (n*SumXY - SumY*SumX)/(n*SumXSquare - SumX**2)
	! Calculate the intercept of the trendline
	intercept = SumY/n - slope*SumX/n
	
END SUBROUTINE CalcSlopeIntercept

REAL FUNCTION CalcCorrCoef(n,xSeries,yValues,slope,intercept)

implicit none

integer, intent(IN) :: n
real, intent(IN), dimension(n) :: xSeries, yValues
real, intent(IN) :: slope, intercept
integer :: counter = 0
real :: AvgValue = 0, error=0, variance=0

	! Calculate the average Temperature for all Temps
	AvgValue = sum(yValues)/real(n)
	
	! Calculate the corrolation coefficient
	error=0
	variance=0
	do counter=1,n
		error = error + (yValues(counter) - (slope*xSeries(counter)+intercept))**2
		variance = variance + (yValues(counter) - AvgValue)**2
	end do

	CalcCorrCoef = 1 - error/variance
	
END FUNCTION CalcCorrCoef