PROGRAM FTCTempDataAnalysis2

! 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) :: AvgMonthlyTemps, TimeIndexes
real, dimension(years) :: YearIndexes, JanTemps, JulTemps
real, dimension(years) :: YearIndexes2, JanTemps2, JulTemps2
integer :: counter = 0, n=years*months
real :: slope = 0, intercept = 0, rsquare=0

    ! 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
	do counter=1,years
		JanTemps(counter)=AvgMonthlyTemps((counter-1)*12+1)
		JulTemps(counter)=AvgMonthlyTemps((counter-1)*12+7)
		YearIndexes(counter)=counter
	end do
	
	! Another way to pull out January's and July's
	JanTemps2 = AvgMonthlyTemps(1:years*months:12)
	JulTemps2 = AvgMonthlyTemps(7:years*months:12)
	YearIndexes2 = floor(TimeIndexes(1:years*months:12)/12.0)+1.0
	do counter=1,years
		print *, 'Year: ', YearIndexes(counter), YearIndexes2(counter), &
		    'Jan Temp: ', JanTemps(counter), JanTemps2(counter), &
			'July Temp: ', JulTemps(counter), JulTemps2(counter)
	end do
	
	! Print statistics for January
	call CalcSlopeIntercept(years, YearIndexes, JanTemps, slope, intercept)
	call CalcCorrCoef(years, YearIndexes, JanTemps, slope, intercept, rsquare)
	print *, 'For January:'
	print *, 'slope = ', slope, ' intercept is ', intercept, ' with r**2= ', rsquare
	
	! Print statistics for July
	call CalcSlopeIntercept(years, YearIndexes, JulTemps, slope, intercept)
	call CalcCorrCoef(years, YearIndexes, JulTemps, slope, intercept, rsquare)
	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

SUBROUTINE CalcCorrCoef(n,xSeries,yValues,slope,intercept,rsquare)

implicit none

integer, intent(IN) :: n
real, intent(IN), dimension(n) :: xSeries, yValues
real, intent(IN) :: slope, intercept
real, intent(OUT) :: rsquare 
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

	rsquare = 1 - error/variance
	
END SUBROUTINE CalcCorrCoef