-1

I have the following script:

import numpy as np
import matplotlib.pyplot as pl

time=np.linspace(0.,1825,1826)
#temp with peak at midpoint of year 
T=23*np.sin((2*np.pi*time)/365)+1
#drainage basin area 
A=8.094e10  
#surface area of pond in mm
SA= 2.023e9 
#hydraulic conductivity in mm/day/mm
K=2e-3
#GW level mm
Hgw=-2000

def precipPond(dayt):
    '''
    Input day with time notation for 365 days... 
    eg time[0] is dayt 1, time[1] is dayt 2, etc.
    '''
    #P varies from 0.2 mm/day to 1.8 mm/day with a peak at 120
    P=np.sin(((2*np.pi*dayt)/365)+120)+0.8
    #15 percent of all P flows into pond (mm)
    infilP=(P*A*0.15)
    return infilP
def evapoTrans(dayT):
    '''
    Input day with time notation for 365 days... 
    eg T[0] is temp at dayT 1, T[1] is temp at dayT 2, etc.
    '''
    #saturated vapor pressure for ET equation
    et=0.611*(np.exp((17.3*dayT)/(dayT+237.3)))
    #evapotranspiration
    ET2=29.8*12*(et)/(dayT+273.2)
    #return factor*ET2
    return ET2

Hp=np.zeros(1826.)
Hp2=[]
precip=[]
evapoT=[]
def waterBudget(HpI,delt):
    '''
    Put in initial water level and change in time.
    '''
    #Save the initial water level in the first position
    Hp[0]=HpI
    steps=(len(time)-1)
    factor=1
    for i in range(steps):
        #Get the mid point value, which will be used to estimate derivative
        HpMid=Hp[i]+(delt/2)*(K*SA*((Hgw-Hp[i])/SA)+precipPond(time[i])+((SA*(precipPond(time[i])/A))*0.85)-(evapoTrans(T[i])*factor*SA))
        #Find the derivative for dh/dt
        Hp[i+1]=Hp[i]+delt*(K*SA*((Hgw-HpMid)/SA)+precipPond(time[i])+((SA*(precipPond(time[i])/A))*0.85)-(evapoTrans(T[i])*factor*SA))
       if any(HpMid < 0):
          factor=0
       else:
          factor=1
    return Hp
ET=evapoTrans(T)
infilP=precipPond(time)
Hp=waterBudget(0,1)

I want the if statement in the loop (within the function waterBudget) to make ET zero when HpMid is below zero. I run this script and it makes factor 1 the whole time even when I've changed the script to make HpMid go below zero (**I didn't include that in here).

Strak
  • 145
  • 1
  • 12
  • 1
    You set `factor`, but it is used nowhere. That's the first problem I see. It also looks like you are trying to call `any` on a float. The rest is too convoluted to understand. Please create a minimum example of what you are trying to do. – timgeb Dec 13 '15 at 02:28

2 Answers2

1

Put:

Hp[i+1]=Hp[i]+delt*( ... )

In the if statement:

HpMid=Hp[i]+(delt/2)*(K*SA*((Hgw-Hp[i])/SA)+precipPond(time[i])+((SA*(precipPond(time[i])/A))*0.85)-(evapoTrans(T[i])*SA))
if HpMid < 0:
    Hp[i+1] = delt*(K*SA*((Hgw-HpMid)/SA)+precipPond(time[i])+((SA*(precipPond(time[i])/A))*0.85)-(evapoTrans(T[i])*SA))
else:
    Hp[i+1]=Hp[i]+delt*(K*SA*((Hgw-HpMid)/SA)+precipPond(time[i])+((SA*(precipPond(time[i])/A))*0.85)-(evapoTrans(T[i])*SA))

you may also want to look at PEP8 for formatting, it will make your code easier to read and thus, easier to debug and reuse.

fiacre
  • 1,150
  • 2
  • 9
  • 26
0

this condition looks wrong if any(HpMid < 0). You can find what any() function does here. In your case HpMid looks like a simple variable holding a double calculated result. so just do if HpMid < 0:

Community
  • 1
  • 1
Xiawei Zhang
  • 1,690
  • 1
  • 11
  • 24
  • If I take out the any(HpMid < 0), I get the following error: ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all() – Strak Dec 13 '15 at 02:34