! basicfun.f90
! module procedure written by Tim Mitchell on 14.12.99
! if you want the slice statistics routines, go look in SliceStats.f90

module BasicFun

implicit none

contains

!*******************************************************************************
! subtract Two from One
! print*, "  > Calc A-B (=1), 100*(A-B)/B (=2), A/B (=3) ?"

function OpTwoFromOne (InputA,InputB,MenuChoice)

real OpTwoFromOne
real, parameter 	:: MissVal = -999.0
real 			:: InputA, InputB
integer 		:: MenuChoice

OpTwoFromOne = MissVal

if (InputA.NE.MissVal.AND.InputB.NE.MissVal) then
     
  if      (MenuChoice.EQ.1) then
     	OpTwoFromOne = InputA - InputB
  else if (MenuChoice.EQ.2) then
     	OpTwoFromOne = 100 * (InputA - InputB) / abs (InputB)
  else if (MenuChoice.EQ.3) then
     	OpTwoFromOne = InputA / InputB
  end if
  
end if

end function OpTwoFromOne

!*******************************************************************************
! calculate weighted mean
! if InputN is unknown it may be set to 'missing'

function OpWeightedMean (Input, Weights, InputN, MissAccept)

real OpWeightedMean
real, pointer, dimension (:) 	:: Input, Weights
real, parameter			:: MissVal = -999.0
real				:: InputSum, WeightSum, TotWeight, MissAccept, Thresh
integer 			:: InputN, XInput

if (InputN.EQ.MissVal) InputN = size (Input,1)

InputSum  = 0.0
WeightSum = 0.0
TotWeight = 0.0

do XInput = 1, InputN
  TotWeight = TotWeight + Weights(XInput)
  
  if (Input(XInput).NE.MissVal.AND.Weights(XInput).NE.MissVal) then
    InputSum  = InputSum  + (Input   (XInput) * Weights (XInput))
    WeightSum = WeightSum +  Weights (XInput)
  end if
end do

Thresh = (100.0 - MissAccept) * TotWeight / 100.0

if (WeightSum.GE.Thresh) then
	OpWeightedMean = InputSum / WeightSum
else
	OpWeightedMean = MissVal
end if

end function OpWeightedMean

!*******************************************************************************
! calculate mean
! if InputN is unknown it may be set to 'missing'

function OpCalcMean (Input, InputN, MissAccept, CallMissVal)

real OpCalcMean
real, pointer, dimension (:) 	:: Input
real, optional			:: CallMissVal
real				:: MissVal
real				:: InputSum, ValidN, MissAccept, Thresh
integer 			:: InputN, XInput

if (present(CallMissVal)) then
  MissVal = CallMissVal
else
  MissVal = -999.0
end if

if (InputN.EQ.MissVal) InputN = size (Input,1)

InputSum = 0.0
ValidN   = 0.0
Thresh   = (100.0 - MissAccept) * InputN / 100.0

do XInput = 1, InputN
  if (Input(XInput).NE.MissVal) then
    InputSum = InputSum + Input(XInput)
    ValidN   = ValidN   + 1.0
  end if
end do

if (ValidN.GE.Thresh) then
	OpCalcMean = InputSum / ValidN
else
	OpCalcMean = MissVal
end if

end function OpCalcMean

!*******************************************************************************
! calculate stdev
! QSampPop: 1=sample-sd, 2=population-sd, 3=sample-variance, 4=pop-variance
! if InputN is unknown it may be set to 'missing'

function OpCalcStDev (Input, InputN, MissAccept, QSampPop, CallMissVal)

real OpCalcStDev
real, pointer, dimension (:) 	:: Input
real, optional			:: CallMissVal
real				:: MissVal
real				:: InputSum, InputSumSq, ValidN, MissAccept, Thresh, Result
integer 			:: InputN, XInput, QSampPop

if (present(CallMissVal)) then
  MissVal = CallMissVal
else
  MissVal = -999.0
end if

if (InputN.EQ.MissVal) InputN = size (Input,1)

InputSum   = 0.0
InputSumSq = 0.0
ValidN     = 0.0
Thresh     = (100.0 - MissAccept) * InputN / 100.0

do XInput = 1, InputN
  if (Input(XInput).NE.MissVal) then
    InputSumSq = InputSumSq + (Input(XInput) ** 2)
    InputSum   = InputSum   +  Input(XInput)
    ValidN     = ValidN     +  1.0
  end if
end do

if (ValidN.GE.Thresh.AND.ValidN.GE.2) then
	Result = (InputSumSq / ValidN) - ( (InputSum / ValidN) ** 2 )
	if (Result.GE.0) then
	  if (QSampPop.EQ.2.OR.QSampPop.EQ.4) Result = Result * ValidN / (ValidN - 1)
	  if (QSampPop.EQ.1.OR.QSampPop.EQ.2) Result = sqrt (Result)
	  OpCalcStDev = Result
	else
	  OpCalcStDev = MissVal
	end if
else
	OpCalcStDev = MissVal
end if

end function OpCalcStDev

!*******************************************************************************
! operate on one value using another
! "  > Divide =1, times =2, minus =3, add =4, sqroot =5, expon =6, abs =7"
! added 30.10.01: 8=log(InputA), 9=e(InputA), 10=ln(InputA)

function OpAwithB (InputA,InputB,MenuChoice)

real OpAwithB
real, parameter 	:: MissVal = -999.0
real 			:: InputA, InputB

integer 		:: MenuChoice, IntegerB

OpAwithB = MissVal

if      (MenuChoice.EQ.1) then
  if (InputA.NE.MissVal.AND.InputB.NE.MissVal.AND.InputB.NE.0) 	OpAwithB = InputA / InputB
else if (MenuChoice.EQ.2) then
  if (InputA.NE.MissVal.AND.InputB.NE.MissVal) 			OpAwithB = InputA * InputB
else if (MenuChoice.EQ.3) then
  if (InputA.NE.MissVal.AND.InputB.NE.MissVal) 			OpAwithB = InputA - InputB
else if (MenuChoice.EQ.4) then
  if (InputA.NE.MissVal.AND.InputB.NE.MissVal) 			OpAwithB = InputA + InputB
else if (MenuChoice.EQ.5) then
  if (InputA.NE.MissVal.AND.InputA.GT.0) 			OpAwithB = sqrt (InputA)
else if (MenuChoice.EQ.6) then
  if (InputA.NE.MissVal.AND.InputB.NE.MissVal) 			OpAwithB = InputA ** InputB
else if (MenuChoice.EQ.7) then
  if (InputA.NE.MissVal) 					OpAwithB = abs (InputA)
else if (MenuChoice.EQ.8) then
  if (InputA.NE.MissVal.AND.InputA.GT.0) 			OpAwithB = log10 (InputA)
else if (MenuChoice.EQ.9) then
  if (InputA.NE.MissVal) 					OpAwithB = exp (InputA)
else if (MenuChoice.EQ.10) then
  if (InputA.NE.MissVal.AND.InputA.GT.0) 			OpAwithB = log (InputA)
end if

end function OpAwithB

!*******************************************************************************

end module BasicFun