text:
# Project Euler problem #112

Problem statement:

Working from left-to-right if no digit is exceeded by the digit to its left it is called an increasing number; for example, 134468.

Similarly if no digit is exceeded by the digit to its right it is called a decreasing number; for example, 66420.

We shall call a positive integer that is neither increasing nor decreasing a "bouncy" number; for example, 155349.

Clearly there cannot be any bouncy numbers below one-hundred, but just over half of the numbers below one-thousand (525) are bouncy. In fact, the least number for which the proportion of bouncy numbers first reaches 50% is 538.

Surprisingly, bouncy numbers become more and more common and by the time we reach 21780 the proportion of bouncy numbers is equal to 90%.

Find the least number for which the proportion of bouncy numbers is exactly 99%.

My code:

import time
start = time.time()

def bouncy(N):
'''Function determines whether arbitrary number N is a bouncy number.'''
number = list(str(N))
n1 = list(str(N))
n2 = list(str(N))
n1.sort()
n2.sort(reverse=True)
if (number != n1 and number != n2):
return True
else:
return False

def find():
start = 1
total = 0
bou = 0
while True:
total += 1
if bouncy(start):
bou += 1
if bou/total >= 0.99:
break
start += 1
return start
print(find())
print(time.time() - start)

What can be improved?

• Just FYI, here is an extended version of the same problem that requires to accept larger inputs (the website provides tests for the extended versions of the first 239 Project Euler problems). For those inputs, a better algorithm is necessary to pass the tests.
– GZ0
Sep 7 '19 at 18:29

# Code structure

The most common and widely accepted structure found in Python scripts from a high-level point of view looks something like:

# imports
import
...

# functions and classes, aka "library code"
class Foo:
...

def bar(batz):
...

if __name__ == "__main__":
# part that is supposed to be run as script
...

Your code on the other hand does not follow this structure since there is some script code at the top (start = time.time()), then come your functions, followed by another block of script code. I highly recommend to follow that structure.

Further reading: Explanation of if __name__ == "__main__": from the official Python documentation.

# The code itself

## bouncy

The most striking aspect of bouncy from a stylistic point is the capital N as parameter name that goes against the official Style Guide for Python Code (often just called PEP8). Usually parameter and variable names are written in snake_case, e.g. small letters, possibly separated by underscores.

With that said, let's look at the actual algorithm:

number = list(str(N))
n1 = list(str(N))
n2 = list(str(N))

Lots of repeated code here. An easy way out would be to do the conversion once, and copy the results afterwards. This can be done using slicing like n1 = number[:], the list constructor n1 = list(number) or the copy module. Also, sorted would create a copy indirectly, so no need to do this manually (see the next point).

The next thing you do in your code is sorting those lists. Once ascending, then descending. Since the second one should by definition be the mirror of the previous one, you can use slicing again to avoid that second sort.

sorted_number = sorted(number)  # <- also does a copy

Oh, and did you know that sorted does work on strings as well and simply returns a new list? (sorted_number = sorted(N))

You could also use the result of the condition directly as return value of your function.

With some other minor changes you might end up with something like

def is_bouncy(number):
'''Function determines whether arbitrary number N is a bouncy number.'''
digits = list(str(number))
sorted_digits = sorted(digits)
return digits != sorted_digits and digits != sorted_digits[::-1]

## find

find is almost as straightforward as it gets. I would advise to change the variable name bou to bouncy or n_bouncy now that the checking function is called is_bouncy. IMHO current would also likely be a better name for start. Also find_bouncy_boundary, maybe even with an optional "expected bounciness percentage" might be worth a thought.

From an optimization point of view, you could also start from 1000 if you're only looking for that 99%, since the task description tells you that there are exactly 525 bouncy numbers below 1000 ;-) But there is likely not much to gain from that.

## The rest

Now all there is left to do is to take the spread out script code, wrap it up in if __name__ == "__main__":, and collect it at the bottom of the script:

if __name__ == "__main__":
start = time.time()
print(find())
print(time.time() - start)

If the script part becomes more complex, this part is often put into a separate main() method.

• number != sorted_digits[::-1] should be digits != sorted_digits[::-1]
– GZ0
Sep 6 '19 at 15:07
• @GZ0: you're correct. Will fix it :-) Sep 6 '19 at 15:28
• Actually, it is better to start from 21780 given from the problem statement that there are 2178 bouncy numbers among them. Sep 7 '19 at 6:43
• @sanyash: That would be 10% then. The task states that exactly 90% are bouncy numbers below 21780, which would be about 19,602 then. Right? Sep 7 '19 at 7:36
• @AlexV yeah, my mistake :) 2178 non-bouncy and 19602 bouncy. Sep 7 '19 at 7:49

First of all, the question asks for the least number which the proportion is exactly 99%. Your code finds the first number that the proportion is more or equal to 99%. That is not quite right. Also, comparing floating-numbers is inaccurate and you should change it to integer comparision for exactness: bou * 100 == total * 99

Secondly, in the bouncy function, list(str(N)) is repeated three times unnecessarily. The twice sort calls are also unnecessary. It can be improved as follows

def bouncy(n):
number = [*str(n)]
sorted_n = sorted(number)
return sorted_n != number and sorted_n[::-1] != number

Thirdly, the overall algorithm is a naive one. If you need a better performance, a better approach is needed, e.g. by exploiting the fact that if $$\n\$$ is a bouncy number, then $$\10n+b\$$ for any $$\b\in[0,9]\$$ is also a bouncy number and does not need to be checked again.

# Speed up the bouncy function

Since you are timing it, I guess speed if of essence. You are sorting a list twice, while you actually don't need to sort it at all. Instead, just check the digits two at a time and see that the difference never change sign. That is O(n) instead of O(n*log(n)) as sorting.

Here's the code. Further changes are highlighted in comments. It only runs slightly faster, 5.8 seconds instead of 8.1 on my machien. My guess is that the reason for the small performance improvement is that Python sorts in C, which makes it comparably fast. But that's just a guess.

import time

# The big change is speeding up this function.
def bouncy(n): # Lower case n
'''Function determines whether arbitrary number N is a bouncy number.'''
# We can index directly into the string. No need to make a list of it.
number = str(n)
allowed_direction = 0
# This could probably be done with cool fancy iterators somehow.
for i in range(0, len(number) - 1):
direction = cmp(number[i], number[i + 1])
if allowed_direction == 0:
allowed_direction = direction
elif direction != 0 and allowed_direction != direction:
return True
return False

def cmp(a, b):
return (a > b) - (a < b)

def find():
# Clearer variable names.
# The start variable was the same as total, so I removed it.
total_count = 100 # Start at 100, no bouncy numbers below anyway!
bouncy_count = 0
while True:
if bouncy(total_count):
bouncy_count += 1
if bouncy_count / total_count >= 0.99:
# No need for a break here. We can just return instead.
total_count += 1

# Cleaner to have this together.
start = time.time()
print(find())
print(time.time() - start)

# You want it faster?

This brute force approach of checking each number is simple, but not fast. When you are at 150 000, you know that the next non bouncy number is going to be 155 555. The 5 555 number in between that you are checking is just a waste of CPU cycles.

Implementing this is significantly harder, and is therefore left as an exercise for the reader.

Edit: Here's a ballpark non solution using the above approach. It runs in 0.15 s on my machine.

# bouncy

You did a good job by extracting the check for bouncyness to a separate function. Plus points for the docstringThe function itself can be a bit better:

You can use the builtin sorted instead of list.sort. And since a string is an iterable too, you don't need the explicit casts to list You can also immediately return the result of the test number != n1 and number != n2

I don't like the names n1 and n2. They are the increasing and decreasing order of digits, so call em that. N

def bouncy(number):
"""Function determines whether arbitrary number N is a bouncy number."""
digits = list(str(number))
increasing = sorted(digits)
decreasing = increasing[::-1]
return digits != decreasing and digits != increasing

I've checked a few variations of this bouncyness check, but this is the fastest I can find.

def pairwise(iterable):
"s -> (s0,s1), (s1,s2), (s2, s3), ..."
a, b = tee(iterable)
next(b, None)
return zip(a, b)

def bouncy_pairwise(number):
"""
tests whether a number is bouncy
"""
digits = str(number)
increasing = decreasing = False
for a, b in pairwise(digits):
if not increasing and a > b:
increasing = True
if not decreasing and a < b:
decreasing = True
if increasing and decreasing:
return True

and

def bouncy_zip(number):
"""
tests whether a number is bouncy
"""
digits = str(number)
increasing = decreasing = False
for a, b in zip(digits, digits[1::]):
if not increasing and a > b:
increasing = True
if not decreasing and a < b:
decreasing = True
if increasing and decreasing:
return True

were faster when the bouncyness was in the beginning of the number, but slower for non-bouncy numbers

# find

• You have the 0.99 hardcoded. I would pass this on as an argument to the list
• Of the double variables total and start, only 1 is needed
• Instead of the break, just return there
• instead of the while True loop, use itertools.count
• no need to shorten bou, you can use bouncy_count or something as a more clear variable name
• You could use the fact that int(True) is 1 and int(False) is 0, to replace if bouncy(i): bouncy_count += 1 by bouncy_count += bouncy(i)

With this result:

def find(fraction = .50):
"""Finds the first number where the ratio of bouncy numbers to all umbers is > fraction"""
bouncy_count = 0
for i in count(1):
bouncy_count += bouncy(i)
if bouncy_count / i > fraction:
return i

# timing

Instead of using time.time, you can use the timeit module

for example:

timeit.repeat("find_bouncy(.9)", globals={"find_bouncy": find_bouncy}, number = 100)
• The checking conditions in bouncy_pairwise can be simplified: if a > b: increasing = True elif a < b: decreasing = True. Another alternative is increasing = increasing or a > b; decreasing = decreasing or a < b.
– GZ0
Sep 6 '19 at 16:13


text:
# p\poly and NP definitions

I have in front of me 1 definition of p\poly and one of NP.

Definition of p\poly:

L E P/poly if there exists a polynomial-time Turing machine M, a polynomial p() and a function h mapping numbers to strings, where |h(n)| <= p(n), such that for all strings x, x E L <=> < x, h(|x|) > is accepted by M.

Definition of NP:

L E NP if there exists a polynomial-time Turing machine M and a polynomial p() such that for all strings x x E L <=> there exists y s.t. |y| <= p(|x|) and < x, y > is accepted by M.

What is the critical difference between the two definitions? For some x1 and x2 with |x1| = |x2| y1 and y2 might be different, but does that mean that p/poly is a subset of NP?

-

There are several "critical differences". The one immediately relevant to your final question is that, in P/poly, $h$ might not even be computable. In the case of NP, the relevant $y$ can be computed, for each $x$ (though it would take exponential time to do so), just by trying each $y$ of the appropriate length for the appropriate number of steps of M. In particular, there are only countably many languages in NP but uncountably many in P/poly.

-

A more illuminating way of thinking about P/poly than the definition you give is in terms of circuit complexity: in fact, P/poly is more like P than like NP (as its name suggests).

It is an exercise to show that for every Turing machine T that computes a function f, there is for each n a circuit (that is, an array of AND, OR and NOT gates) that takes all inputs x of length n and computes f(x). Moreover, the size of the circuit differs from the time taken by T by at most a polynomial factor.

Thus, any Turing machine can be replaced by a sequence of circuits, one for each input length. Moreover, this sequence of circuits can itself be generated by an algorithm.

If you take an arbitrary sequence of circuits (as opposed to one generated by an algorithm) taking increasing input sizes, then the resulting function is said to belong to P/poly. That is, a function in P/poly is like a function in P but without the property of uniformity -- the property that the circuits are all generated by the same algorithm. It is for that reason that you can get non-computable functions, that there are uncountably many of them, etc. etc.

The additional string in your formal definition can be thought of as encoding a circuit.

If you think about things this way it is much more obvious that the definition is not the same as that of NP.

-

The definitions that you give are very formal and tend to hide the conceptual interpretation of these complexity classes.

Imagine putting yourself in the position of the computer or Turing machine that is allowed to work for a polynomial length of time. (It's a reasonable model if you believe that people are evolved computers.) The string h in the definition of P/poly is called an advice string. You can say that it is provided by an angel who wants to help you produce the correct answer. If h could depend on the input x generally, then the angel could simply tell you the answer. So instead h only depends on the length of the input |x|.

Meanwhile in the class NP, the string y is called a proof or a certificate. It is provided by a Merlin who wants you to answer "yes" to the decision problem, whether or not it is the correct answer. It can depend on the specific input x and not just on its length |x|.

As Andreas says, a problem in P/poly does not even need to be computable. But I can mention a realistic class of problems that are in P/poly and are computable. The class BPP, which is anything computable probably correctly in polynomial time with a random number generator, is in P/poly, by a derandomization theorem. The theorem says that the angel can provide pseudorandom numbers that, together with repeated trials, make the randomized algorithm reliable as a non-randomized (but advice-dependent) algorithm.

On the the other hand, an NP-complete problem such as the Hamiltonian circuit is easily seen to fit the Merlin model. If Merlin wants to convince you that there is a Hamiltonian circuit, he only needs to show you what it is. If his certificate y isn't a Hamiltonian circuit, then you should be skeptical and answer "no".

-

format: plain


text:
histogram¶

astropy.stats.histogram(a, bins=10, range=None, weights=None, **kwargs)[source]

Enhanced histogram function, providing adaptive binnings

This is a histogram function that enables the use of more sophisticated algorithms for determining bins. Aside from the bins argument allowing a string specified how bins are computed, the parameters are the same as numpy.histogram().

Parameters
aarray_like

array of data to be histogrammed

binsint, list, or str, optional

If bins is a string, then it must be one of:

• ‘blocks’ : use bayesian blocks for dynamic bin widths

• ‘knuth’ : use Knuth’s rule to determine bins

• ‘scott’ : use Scott’s rule to determine bins

• ‘freedman’ : use the Freedman-Diaconis rule to determine bins

rangetuple or None, optional

the minimum and maximum range for the histogram. If not specified, it will be (x.min(), x.max())

weightsarray_like, optional

An array the same shape as a. If given, the histogram accumulates the value of the weight corresponding to a instead of returning the count of values. This argument does not affect determination of bin edges.

other keyword arguments are described in numpy.histogram().
Returns
histarray

The values of the histogram. See density and weights for a description of the possible semantics.

bin_edgesarray of dtype float

Return the bin edges (length(hist)+1).


text:
No such thing as a small change PerlMonks

### Getting Cookies

by PerlSufi (Friar)
 on Jul 05, 2013 at 15:26 UTC ( #1042736=perlquestion: print w/ replies, xml ) Need Help??
PerlSufi has asked for the wisdom of the Perl Monks concerning the following question:

Hello monks,
I want to extract cookies from a website. Particularly, JSESSIONID cookie. However, when I run the code below, nothing turns up in my file except '#LWP-Cookies-1.0'.
use LWP::UserAgent;
use HTTP::Cookies;

my $cookie_jar = HTTP::Cookies->new( file => 'C:\path\to\cookies.dat', + autosave=> 1); my$browser = LWP::UserAgent->new( );

$browser->cookie_jar($cookie_jar );

$browser->get('https://example.com'); [download] Any insight into what I'm doing wrong is greatly appreciated! Replies are listed 'Best First'. Re: Getting Cookies by kennethk (Abbot) on Jul 05, 2013 at 16:02 UTC Your code above is correct, assuming that your cookie isn't timing out (i.e. isn't a session cookie). Under normal circumstances, session cookies are discarded immediately upon exit from your 'browser'. If you want to force holding onto the cookie, use ignore_discard in your constructor: my$cookie_jar = HTTP::Cookies->new(file => 'C:\path\to\cookies.dat',
autosave=> 1,
ignore_discard => 1);
[download]
The flag is documented in HTTP::Cookies, and read HTTP_cookie#Session_cookie to learn more about the behavior.

#11929 First ask yourself How would I do this without a computer?' Then have the computer do it the same way.

Thanks, kennethk. That worked for me. :)
The cookie is different from what my browser gives me so I'm working on how to get that now
Firebug, the free Firefox plug-in, makes that easy.

I'd also mention that, instead of storing cookies to disk, you can use Data::Dumper to look at your cookie jar contents.

#11929 First ask yourself How would I do this without a computer?' Then have the computer do it the same way.

Re: Getting Cookies (webbug)
by Anonymous Monk on Jul 06, 2013 at 06:46 UTC

Well, for some reason you think https://example.com/ will set some cookies for you, and not, for example, an image embedded inside html returned by https://example.com/ or javascript embedded inside or css embedded inside or one or more combinations of both

Also, could be https://example.com/ doesn't like the headers you're sending, so it doesn't give you a cookie

Also, it could be a proxy of some sort between you and https://example.com/ interfering

Also, it could be some firewall or antivirus or somesuch interfering on your machine (or on proxy between you and https://example.com/ )

Also, it could be anything :)

Log In?
 Username: Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: perlquestion [id://1042736]
Front-paged by Arunbear
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others studying the Monastery: (5)
As of 2016-07-24 06:08 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
What is your favorite alternate name for a (specific) keyboard key?

Results (221 votes). Check out past polls.


text:
### Top 9 Arxiv Papers Today in Analysis Of Pdes

##### #1. Justification of the Point Vortex Approximation for Modified Surface Quasi-Geostrophic Equations
###### Matthew Rosenzweig
In this paper, we give a rigorous justification of the point vortex approximation to the family of modified surface quasi-geostrophic (mSQG) equations globally in time in both the inviscid and vanishing dissipative cases. This result completes the justification for the remaining range of the mSQG family unaddressed by Geldhauser and Romito in the article arXiv:1812.05166, in the case of identically signed vortices.
more | pdf | html
None.
###### Tweets
mathAPb: Matthew Rosenzweig : Justification of the Point Vortex Approximation for Modified Surface Quasi-Geostrophic Equations https://t.co/UA2ExZSpbI https://t.co/aqjGF8yvhh
MathPHYPapers: Justification of the Point Vortex Approximation for Modified Surface Quasi-Geostrophic Equations. https://t.co/zWK7GzbVV1
None.
None.
###### Other stats
Sample Sizes : None.
Authors: 1
Total Words: 16779
Unqiue Words: 3079

##### #2. The Gabor wave front set of compactly supported distributions
###### Patrik Wahlberg
We show that the Gabor wave front set of a compactly supported distribution equals zero times the projection on the second variable of the classical wave front set.
more | pdf | html
None.
###### Tweets
mathAPb: Patrik Wahlberg : The Gabor wave front set of compactly supported distributions https://t.co/RJDPgLdzE9 https://t.co/KqVgTd1VfL
None.
None.
###### Other stats
Sample Sizes : None.
Authors: 1
Total Words: 4468
Unqiue Words: 1166

##### #3. On the Well-posedness of a Nonlinear Fourth-Order Extension of Richards' Equation
###### Alaa Armiti-Juber, Christian Rohde
We study a nonlinear fourth-order extension of Richards' equation that describes infiltration processes in unsaturated soils. We prove the well-posedness of the fourth-order equation by first applying Kirchhoff's transformation to linearize the higher-order terms. The transformed equation is then discretized in time and space and a set of a priori estimates is established. These allow, by means of compactness theorems, extracting a unique weak solution. Finally, we use the inverse of Kirchhoff's transformation to prove the well-posedness of the original equation.
more | pdf | html
None.
###### Tweets
mathAPb: Alaa Armiti-Juber, Christian Rohde : On the Well-posedness of a Nonlinear Fourth-Order Extension of Richards' Equation https://t.co/PcCaYGTP02 https://t.co/4NnSbYmwSk
None.
None.
###### Other stats
Sample Sizes : None.
Authors: 2
Total Words: 0
Unqiue Words: 0

##### #4. Asymptotic flatness of Morrey extremals
###### Ryan Hynd, Francis Seuffert
We study the limiting behavior as $|x|\rightarrow \infty$ of extremal functions $u$ for Morrey's inequality on $\mathbb{R}^n$. In particular, we compute the limit of $u(x)$ as $|x|\rightarrow \infty$ and show $|x||Du(x)|$ tends to $0$. To this end, we exploit the fact that extremals are uniformly bounded and that they each satisfy a PDE of the form $-\Delta_pu=c(\delta_{x_0}-\delta_{y_0})$ for some $c\in\mathbb{R}$ and distinct $x_0,y_0\in\mathbb{R}^n$. More generally, we explain how to deduce the asymptotic flatness of bounded $p$-harmonic functions on exterior domains of $\mathbb{R}^n$ for $p>n$.
more | pdf | html
###### Tweets
mathAPb: Ryan Hynd, Francis Seuffert : Asymptotic flatness of Morrey extremals https://t.co/gy634GNZAJ https://t.co/C8UJ2DIcY5
None.
None.
###### Other stats
Sample Sizes : None.
Authors: 2
Total Words: 5954
Unqiue Words: 1375

##### #5. An example of failure of stochastic homogenization for viscous Hamilton-Jacobi equations without convexity
###### William M. Feldman, Jean-Baptiste Fermanian, Bruno Ziliotto
We give an example of the failure of homogenization for a viscous Hamilton-Jacobi equation with non-convex Hamiltonian.
more | pdf | html
None.
###### Tweets
mathAPb: William M. Feldman, Jean-Baptiste Fermanian, Bruno Ziliotto : An example of failure of stochastic homogenization for viscous Hamilton-Jacobi equations without convexity https://t.co/QOcl1o8vKo https://t.co/sg19icH5Xi
None.
None.
###### Other stats
Sample Sizes : None.
Authors: 3
Total Words: 0
Unqiue Words: 0

##### #6. Kramers-Fokker-Planck operators with homogeneous potentials
###### Mona Ben Said
In this article we establish a global subelliptic estimate for Kramers-Fokker-Planck operators with homogeneous potentials $V(q)$ under some conditions, involving in particular the control of the eigenvalues of the Hessian matrix of the potential. Namely, this work presents a different approach from the one in [Ben], in which the case $V (q_1, q _2) =-q ^2_1 (q^ 2_1+q^2_2) ^n$ was already treated only for $n=1.$ With this article, after the former one dealing with non homogeneous polynomial potentials, we conclude the analysis of all the examples of degenerate ellipticity at infinty presented in the framework of Witten Laplacian by Helffer and Nier in [HeNi]. Like in [Ben], our subelliptic lower bounds are the optimal ones up to some logarithmic correction.
more | pdf | html
None.
###### Tweets
mathAPb: Mona Ben Said (LAGA) : Kramers-Fokker-Planck operators with homogeneous potentials https://t.co/J0zTAPAH0E https://t.co/IysT8aKxFP
None.
None.
###### Other stats
Sample Sizes : None.
Authors: 1
Total Words: 6239
Unqiue Words: 1246

##### #7. Representation of solutions to wave equations with profile functions
###### Agnes Lamacz, Ben Schweizer
Solutions to the wave equation with constant coefficients in $\mathbb{R}^d$ can be represented explicitly in Fourier space. We investigate a reconstruction formula, which provides an approximation of solutions $u(.,t)$ to initial data $u_0(.)$ for large times. The reconstruction consists of three steps: 1) Given $u_0$, initial data for a profile equation are extracted. 2) A profile evolution equation determines the shape of the profile at time $\tau = \varepsilon^2 t$. 3) A shell reconstruction operator transforms the profile to a function on $\mathbb{R}^d$. The sketched construction simplifies the wave equation, since only a one-dimensional problem in an $O(1)$ time span has to be solved. We prove that the construction provides a good approximation to the wave evolution operator for times $t$ of order $\varepsilon^{-2}$.
more | pdf | html
###### Tweets
mathAPb: Agnes Lamacz, Ben Schweizer : Representation of solutions to wave equations with profile functions https://t.co/ztQu5jwaby https://t.co/ms4fBgj6SR
None.
None.
###### Other stats
Sample Sizes : None.
Authors: 2
Total Words: 8193
Unqiue Words: 1573

##### #8. $W^{s,\frac{n}{s}}$-maps with positive distributional Jacobians
###### Siran Li, Armin Schikorra
We extend the well-known result that any $f \in W^{1,n}(\Omega,\mathbb{R}^n)$, $\Omega \subset \mathbb{R}^n$ with strictly positive Jacobian is actually continuous: it is also true for fractional Sobolev spaces $W^{s,\frac{n}{s}}(\Omega)$ for any $s \geq \frac{n}{n+1}$, where the sign condition on the Jacobian is understood in a distributional sense. Along the way we also obtain extensions to fractional Sobolev spaces $W^{s,\frac{n}{s}}$ of the degree estimates known for $W^{1,n}$-maps with positive or non-negative Jacobian, such as the sense-preserving property.
more | pdf | html
None.
###### Tweets
mathAPb: Siran Li, Armin Schikorra : $W^{s,\frac{n}{s}}$-maps with positive distributional Jacobians https://t.co/ju5hXcBFYK https://t.co/sNBFKa6S6O
None.
None.
###### Other stats
Sample Sizes : None.
Authors: 2
Total Words: 0
Unqiue Words: 0

##### #9. Large time asymptotics for a cubic nonlinear Schrödinger system in one space dimension
###### Chunhua Li, Yoshinori Nishii, Yuji Sagawa, Hideaki Sunagawa
We consider a two-component system of cubic nonlinear Schr\"odinger equations in one space dimension. We show that each component of the solutions to this system behaves like a free solution in the large time, but there is a strong restriction between the profiles of them. This turns out to be a consequence of non-trivial long-range nonlinear interactions.
more | pdf | html
None.
###### Tweets
mathAPb: Chunhua Li, Yoshinori Nishii, Yuji Sagawa, Hideaki Sunagawa : Large time asymptotics for a cubic nonlinear Schrödinger system in one space dimension https://t.co/Fkvxy72FsN https://t.co/2CRQGjR8fj
None.
None.
###### Other stats
Sample Sizes : None.
Authors: 4
Total Words: 5065
Unqiue Words: 1368

Assert is a website where the best academic papers on arXiv (computer science, math, physics), bioRxiv (biology), BITSS (reproducibility), EarthArXiv (earth science), engrXiv (engineering), LawArXiv (law), PsyArXiv (psychology), SocArXiv (social science), and SportRxiv (sport research) bubble to the top each day.

Papers are scored (in real-time) based on how verifiable they are (as determined by their Github repos) and how interesting they are (based on Twitter).

To see top papers, follow us on twitter @assertpub_ (arXiv), @assert_pub (bioRxiv), and @assertpub_dev (everything else).

To see beautiful figures extracted from papers, follow us on Instagram.

Tracking 128,326 papers.

###### Search
Sort results based on if they are interesting or reproducible.
Interesting
Reproducible
Online
###### Stats
Tracking 128,326 papers.


text:
# What is the distance between (0,0) and (-15,36)?

May 13, 2016

$39$

#### Explanation:

From Pythagoras theorem, we get the following formula for the distance between points $\left({x}_{1} , {y}_{1}\right)$ and $\left({x}_{2} , {y}_{2}\right)$ in the plane:

$d = \sqrt{{\left({x}_{2} - {x}_{1}\right)}^{2} + {\left({y}_{2} - {y}_{1}\right)}^{2}}$

In our example, $\left({x}_{1} , {y}_{1}\right) = \left(0 , 0\right)$ and $\left({x}_{2} , {y}_{2}\right) = \left(- 15 , 36\right)$, giving us:

$d = \sqrt{{\left(- 15 - 0\right)}^{2} + {\left(36 - 0\right)}^{2}}$

$= \sqrt{{\left(- 15\right)}^{2} + {36}^{2}}$

$= \sqrt{225 + 1296}$

$= \sqrt{1521} = 39$


text:
# help working out if 2 spheres overlap / colided

This topic is 4761 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

i have been looking on gamasutra,at the sphere on sphere sweep tests, but i cant figure out the code, partly as its in c++ ( i only know c#) and 2nd that im not fantastic at maths. can anyone explain better how i can work out if spheres have colided?

##### Share on other sites
Moved to "Math and physics".

##### Share on other sites
Quote:
 can anyone explain better how i can work out if spheres have colided?

Maybe I'm over simplifying it, but given a center and radius for each sphere, the following logic is all thats needed:

1. Calculate the distance between the two center points:
dist = square_root( a*a + b*b + c*c )
where:
- a = (sphere1_x - sphere2_x)
- b = (sphere1_y - sphere2_y)
- c = (sphere1_z - sphere2_z)

2. if the distance is GREATER THAN the sum of the two radii then the spheres are not touching/intersecting

3. if the distance is EQUAL TO the sum of the two radii they are just touching on the surfaces.

4. if the distance is LESS THEN the sum of the two radii they are intersecting/overlapping each other.

I'm sure you can convert that into relatively simple C# code [grin]

hth
Jack

##### Share on other sites

What part of the code/algorithm are you having trouble with? The algorithm is fairly straightforward and involves solving a quadratic equation to find the first time, if any, of contact between the spheres. If you can tell us what you're uncertain about, someone here will be able to clarify things for you.

##### Share on other sites
Quote:

oops! my bad... thought it sounded a bit too simple - I usually don't get my hands dirty in M&P for this reason [grin]..

Jack

##### Share on other sites
It is easy enough to see that sweep sphere test is the same as ray-sphere intersection. Let one sphere is stationary, other is moving along some ray. Let radiuses is r1,r2 respectively. Then, spheres collide when distance between centers <r1+r2 . Same as intersection of ray with sphere that have radius r1+r2.

If both spheres is moving, you need to do everything relatively to one of spheres.

I have posted my ray-sphere intersection code to this forum about 7 times.
//// Ray-sphere intersection.// p=(ray origin position - sphere position),// d=ray direction,// r=sphere radius,// Output:// i1=first intersection distance,// i2=second intersection distance// i1<=i2// i1>=0// returns true if intersection found,false otherwise.//bool RaySphereIntersect(const Vec3 &p, const Vec3 &d,double r, double &i1, double &i2){        double det,b;        b = -DotProduct(p,d);        det=(b*b) - DotProduct(p,p) + r*r;        if (det<0){                return false;        }        det= sqrt(det);        i1= b - det;        i2= b + det;        // intersecting with ray?        if(i2<0) return false;        if(i1<0)i1=0;        return true;};

note that direction must be normalized.
If you call this with
// doing everything relatively to first sphere
Vec3 V=Velocity2-Velocity1;
double l=Length(V);
V*=1/l;
double i1,i2
i1/=l;
you'll get i1=time of first collision. i1=0 if already collided.

You can optimize code if you like. I think, that
// // Warning: Non-tested function.//// Ray-sphere intersection that works as moving point-sphere collision.// p=(point start position - sphere center position),// v=point velocity,// r=sphere radius,// Output:// i1=first intersection time,// i2=second intersection time// i1<=i2// i1>=0// returns true if intersection found,false otherwise.//bool RaySphereIntersectTime(const Vec3 &p, const Vec3 &v,double r, double &i1, double &i2){        double det,b,k;        k = DotProduct(v,v);          b = -DotProduct(p,v)/k;                det=b*b - (DotProduct(p,p) - r*r)/k;        if (det<0){                return false;        }        det= sqrt(det);        i1= b - det;        i2= b + det;        // intersecting with ray?        if(i2<0) return false;        if(i1<0)i1=0;        return true;};


text:
# 7th Class Mathematics Integers Properties of Integers

Properties of Integers

Category : 7th Class

### Properties of Integers

There are four properties of Integers

• Properties of Subtraction
• Properties of Multiplication
• Properties of Division

Closure Property

The sum of two integers is always an integer.

For any two integer's p and q, $\Rightarrow p+q$ is an integer.

For example $\left( -\text{ }4 \right)+3=\left( -\text{ }1 \right),$ Which is an integer.

Commutativity

Sum of two integers is commutative.

For any two integers p and q,$~\Rightarrow p+q=q+p,$ for example $7+\left( -6 \right)=1$and $\left( -\text{ }6 \right)+7=1,$ therefore $7\text{ }+\text{ }\left( -\text{ }6 \right)\text{ }=\text{ }\left( -6 \right)+7$

Associativity

We addition of integers is associative. For any three integers p, q and r, $p+\left( q+r \right)=\left( p+q \right)+r,$

$\left( -4 \right)+\left\{ \left( -8 \right)+6 \right\}=\left\{ \left( -4 \right)+\left( -8 \right) \right\}+6$

Solution:

L.H.S $=\left( -4 \right)+\left\{ \left( -\text{ }8 \right)+6 \right\}R.H.S=\left\{ \left( -4 \right)+\left( -8 \right) \right\}+6$

$=(-4)+(-2)=(-12)+6=(-6)=(-6)$

Additive identity is the number which when added to an integer gives the same integer. In case of addition "0" is the identity for integers. For example:

(i)  $~\left( -8 \right)+0=\left( -8 \right)$

(ii)   $~0+\left( -37 \right)=\left( -37 \right)$

(iii)   $31+0=37$

The additive inverse of an integer is that integer which when added to the given integer gives the additive identity $0.\left( -p \right)$ is the additive inverse of p, i.e. $p+\left( -p \right)$ = 0. We get the additive inverse by taking the negative of the given integer. For example ( - 10) is additive inverse of 10 and 20 is the additive inverse of (-20)

Properties of Subtraction

Closure Property

Integers are closed under subtraction. If a and b are two integers then a - b is also an integer.

For example

(i) $7-9=-2$ is an integer

(ii) $\left( -8 \right)-\left( -15 \right)=7$ is an integer

Subtraction is neither Commutative nor Associative

For any two numbers a and $b,a-b\ne b-a$ and for any three integers a, b, c we have $a-(b-c)\ne (a-b)-c$

(i)  $9-3\ne 3-9$

(ii) $9-(3-2)\ne (9-3)-2$

A bird is flying at the height of 750 m above the lake surface. At a particular point it is exactly above a fish swimming 120 m below the lake surface. What is the vertical separation between bird and fish?

(a) 150m

(b) 350m

(c) 870m

(d) 630m

(e) Data is insufficient

Explanation

Distance between bird and fish | 750 m|+|-120 m| = 750 m+120 m = 870 m

Therefore, option (C) is correct and rest of the options is incorrect.

In London temperature was $-\text{ }9{}^\circ c$on Monday and then dropped by $2{}^\circ c$ on next day. On Wednesday it raised by again. The temperature on Wednesday was:

(a)$7{}^\circ c$

(b) $-4{}^\circ c$

(c)$-11{}^\circ c$

(d) 18°c

(e) None of these

John goes 17 km towards north from a point P to another point Q. From Q he moves 27 km towards south along the same road. If the distance towards north is represented by positive integer then which integer will represent his final distance from P?

(a) $-10km$

(b) 44km

(c)$~-34km~$

(d) 10km

(e) None of these

There are 40 questions in the test booklet of Mathematics Olympiad. If Stephen scores 80 marks when 5 marks are awarded for each correct answer and 3 marks is deducted for each incorrect answer then find the number of correct answer attempted by Stephen.

(a) 16

(b) 20

(c) 28

(d) 25

(e) None of these

If X and Y are two positive integers so that X > Y and P, Q are the additive inverse of X and Y respectively, then which of the following is correct?

(a)$~X+Y+P-Q=Y~$

(b)$~\left( X+Y \right)-\left( P+Q \right)=0$

(c)$2\left( X+Y \right)+\left( P+Q \right)=-\left( P+Q \right)$

(d) $\left( X+Y \right)+\left( P+Q \right)=2\left( P+Q \right)$

(e) None of these

Properties of Multiplication

Closure Property

For any two integers 'm' and 'n', the product ?m' and 'n' is always an integer.

For example:

(i)  $8\times 3=24\text{ },$which is an integer.

(ii)  $~7\times \left( -5 \right)=-35,$which is an integer.

Commutativity

Multiplication of integers is commutative.

For any two integers m and n, $m\times n=n\times m$

For example: $8x\left( -\text{ }7 \right)=-\left( 8\text{ }x\text{ }7 \right)=-56$  and $\left( -\text{ }7 \right)\times 8=-\left( 7\text{ }\times \text{ }8 \right)=-56$

Associativity

The multiplication of integers is associative.

For any three integers m, n and $p,\times m\times \left( n\times p \right)=\left( m\times n \right)\times \text{ }p,$for example $\left( -12 \right)\times \left\{ 5\text{ }x\left( -2 \right) \right\}=\left( -12 \right)\times \left( -10 \right)=120$ and $\left\{ \left( -12 \right)\times 5 \right\}\times \left( -2 \right)=\left( -60 \right)\times \left( -2 \right)=120.$ $\therefore \left( -\text{ }12 \right)\times \left\{ 5\times \left( -2 \right) \right\}=\left\{ \left( -12 \right)\times 5 \right\}\times \left( -2 \right)$

Multiplicative Identity

The multiplicative identity of a number is that number which when multiplied with the given integer results the same integers. "1" is the multiplicative identity for any integers. For any integers 'm' we can write $1\times m=m=m\times 1$

Multiplicative Inverse

An inverse element usually does not exist at all. For example "0" is an integer whose inverse does not exist.

Distributive Property of Multiplication on Addition

For any three integers m, n, p, we have:

(i)    $m\times \left( n+p \right)=\left( m\times n \right)+\left( m\times p \right)$

(ii)  $\left( n+p \right)\times m=n\times m+n\times p$

Properties of Division

For any two integer p and q

• For any two positive integer's p and q, p - q is not necessarily an integer. $25\div 7,\left( -\text{ }125 \right)\div 23$
• $p\div p=1,$if p is any integer other than 0.
• For any integer p, we have $p\div 1=p$
• There are three integers p, q and r other than "0" then

(i)  $p<q\Rightarrow p\div r<q\div r,$if r is positive.

(ii)  $p<q\Rightarrow p\div r>q\div r,$if r is negative.

What feast value must be given to '*' so that the number 97215*6 is divisible by 11 ?

(a) 11

(b) 12

(c) 3

(d) 15

(e) None of these

Explanation

$\left( 9+2+5+6 \right)-\left( 7+1+x \right)=14-x,$must be divisible by 11. So, $x=3.$

What least value must be given to '*' so that the number 91876*2 is divisible by 8?

(A) 1

(b) 2

(c) 3

(d) 4

(e) None of these

Explanation

By hit and trial, we find that 9187632 is divisible by 8. So '*' must be replaced by 3.

Which one of the following numbers is exactly divisible by 99?

(a) 3572404

(b) 135792

(c) 913464

(d) 114345

(e) None of these

Explanation

We know that if a number is divisible by 9 and 11 both then that number must be divisible by 99. Clearly 114345 is divisible by 9 and 11. So, it is divisible by 99.

$\frac{392}{\sqrt{?}}=28$

(a) 144

(b) 196

(c) 24

(d) 48

(e) None of these

Explanation

let, $\frac{392}{\sqrt{x}}=28$

Therefore, $28\sqrt{x}=392${By cross multiplication}

$\sqrt{x}=\frac{392}{28}\Rightarrow x={{(14)}^{2}}=196$

The largest number of four digits exactly divisible by 88 is:

(a) 9768

(b) 8888

(c) 9988

(d) 9944

(e) None of these

Explanation

Largest number of four digits = 9999.

Required number $=\text{ }9999\text{ }-$ {Remainder of $\left( 9999\div 88 \right)\}$

$\therefore$Required number $=9999-55=9944$

The least number of five digit exactly divisible by 456 is?

(a) 10142

(b) 10232

(c) 10032

(d) 10012

(e) None of these

Explanation

Least number of five digits = 10000

Required number $=10000+(456-$remainder of division, 10000 by 456) $=\text{ }10000\text{ }+\text{ }\left( 456-424 \right)=10000+32=10032.$

What least number must be subtracted from 13601 to get a number exactly divisible by 87?

(a) 49

(b) 23

(c) 29

(d) 31

(e) None of these

What least number must be added to 1056 to get a number exactly divisible by 23?

(a) 21

(b) 25

(c) 13

(d) 2

(e) None of these

What least value must be given to '*' so that the number 6135*2 is exactly divisible by 9?

(a) 0

(b) 1

(c) 12

(d) 13

(e) None of these


text:
# [Grammar]Which one is right? (easy)

Status
Not open for further replies.

#### DSoldier

##### Member
Is it your first time here Alan?

or

Is this your first time here Alan?

Which one of them is the right one? :shock:

Thank you all in advance! :up:

*ilovethisforum*

#### Barb_D

##### Moderator
Staff member
Hi, and welcome to Using English.

Either one, although I think I'd be more likely to use "this" myself.

You should put a comma before "Alan."

#### sarat_106

##### Key Member
Hi, and welcome to Using English.

Either one, although I think I'd be more likely to use "this" myself.

You should put a comma before "Alan."

Yes, I agree. Both are potentially correct. The only difference is that ‘this’ refer to a thing in general, much more directly than "it" without a previous reference. So, use ‘it’ with a prior reference; as:
Welcome to this nice place. Is it your first time here, Alan?

Status
Not open for further replies.


text:
# Implementation of the Thornthwaite-Mather procedure to map groundwater recharge

Author(s): guiattard

Groundwater recharge represents the amount of water coming from precipitation reaching the groundwater table. Its determination helps to better understand the available/renewable groundwater in watersheds and the shape of groundwater flow systems.

One of the simplest methods to estimate groundwater recharge is the Thornthwaite-Mather procedure (Steenhuis and Van Der Molen, 1986). This procedure was published by Thornthwaite and Mather (1955, 1957). The idea of this procedure is to calculate the water balance in the root zone of the soil where water can be (1) evaporated into the atmosphere under the effect of heat, (2) transpired by vegetation, (3) stored by the soil, and eventually (4) infiltrated when stored water exceeds the field capacity.

This procedures relies on several parameters and variables described as follows:

• information about soil texture (e.g. sand and clay content) to describe the hydraulic properties of the soil and its capacity to store/infiltrate,
• meteorological records: precipitation and potential evapotranspiration.

Of course groundwater recharge can be influenced by many other factors such as the slope of the terrain, the snow cover, the variability of the crop/land cover and the irrigation. In the following these aspects are not taken into account.

In the first part of the tutorial, the Earth Engine python API will be initialized, some useful libraries will be imported, and the location/period of interest will be defined.

In the second part, OpenLandMap datasets related to soil properties will be explored. The wilting point and field capacity of the soil will be calculated by applying some mathematical expressions to multiple images.

In the third part, evapotranspiration and precipitation datasets will be imported. A function will be defined to resample the time resolution of an ee.ImageCollection and to homogenize time index of both datasets. Both datasets will then be combined into one.

In the fourth and final part, the Thornthwaite-Mather(TM) procedure will be implemented by iterating over the meteorological ee.ImageCollection. Finally, a comparison between groundwater recharge in two places will be described and the resulting mean annual groundwater recharge will be displayed over France.

### Run me first

#### Earth Engine API

First of all, run the following cell to initialize the API. The output will contain instructions on how to grant this notebook access to Earth Engine using your account.

import ee

# Trigger the authentication flow.
ee.Authenticate()

# Initialize the library.
ee.Initialize()

#### Other libraries

Import other libraries/modules used in this notebook.

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import folium
import pprint
import branca.colormap as cm

#### Some input parameters

We additionally define some parameters used to evaluate the results of our implementation. Particularly:

• the period of interest to get meteorological records,
• a location of interest based on longitude and latitude coordinates. In the following, the point of interest is located in a productive agricultural region which is about 30 kilometers outside of the city of Lyon (France). This point is used to evaluate and illustrate the progress of the described procedure.
# Initial date of interest (inclusive).
i_date = "2015-01-01"

# Final date of interest (exclusive).
f_date = "2020-01-01"

# Define the location of interest with a point.
lon = 5.145041
lat = 45.772439
poi = ee.Geometry.Point(lon, lat)

# A nominal scale in meters of the projection to work in [in meters].
scale = 1000

### From soil texture to hydraulic properties

#### Definitions

Two hydraulic properties of soil are commonly used in the TM procedure:

• the wilting point represents the point below what water cannot be extracted by plant roots,
• the field capacity represents the point after which water cannot be stored by soil any more. After that point, gravitational forces become too high and water starts to infiltrate the lower levels.

Some equations given by Saxton & Rawls (2006) are used to link both parameters to the texture of the soil. The calculation of water content at wilting point $$θ_{WP}$$ can be done as follows:

$$\theta_{WP}= \theta_{1500t} + (0.14 \theta_{1500t} - 0.002)$$

with:

$$\theta_{1500t} = -0.024 S + 0.487 C + 0.006 OM + 0.005(S \times OM) - 0.013 (C \times OM) + 0.068 (S \times C) + 0.031$$

where:

• $$S$$: represents the sand content of the soil (mass percentage),
• $$C$$: represents the clay content of the soil (mass percentage),
• $$OM$$: represents the organic matter content of the soil (mass percentage).

Similarly, the calculation of the water content at field capacity $$θ_{FC}$$ can be done as follows:

$$\theta_{FC} = \theta_{33t} + (1.283 \theta_{33t}^{2} - 0.374 \theta_{33t}-0.15)$$

with:

$$\theta_{33t} = -0.251 S + 0.195 C + 0.011 OM + 0.006 (S \times OM) - 0.027 (C \times OM) + 0.452 (S \times C) + 0.299$$

#### Determination of soil texture and properties

In the following, OpenLandMap datasets are used to describe clay, sand and organic carbon content of soil. A global dataset of soil water content at the field capacity with a resolution of 250 m has been made available by Hengl & Gupta (2019). However, up to now, there is no dataset dedicated to the water content of soil at the wilting point. Consequently, in the following, both parameters will be determined considering the previous equations and using the global datasets giving the sand, clay and organic matter contents of the soil. According to the description, these datasets are based on machine learning predictions from global compilation of soil profiles and samples. Processing steps are described in detail here. The information (clay, sand content, etc.) is given at 6 standard depths (0, 10, 30, 60, 100 and 200 cm) at 250 m resolution.

These standard depths and associated bands are defined into a list as follows:

# Soil depths [in cm] where we have data.
olm_depths = [0, 10, 30, 60, 100, 200]

# Names of bands associated with reference depths.
olm_bands = ["b" + str(sd) for sd in olm_depths]

We now define a function to get the ee.Image associated to the parameter we are interested in (e.g. sand, clay, organic carbon content, etc.).

def get_soil_prop(param):
"""
This function returns soil properties image
param (str): must be one of:
"sand"     - Sand fraction
"clay"     - Clay fraction
"orgc"     - Organic Carbon fraction
"""
if param == "sand":  # Sand fraction [%w]
snippet = "OpenLandMap/SOL/SOL_SAND-WFRACTION_USDA-3A1A1A_M/v02"
# Define the scale factor in accordance with the dataset description.
scale_factor = 1 * 0.01

elif param == "clay":  # Clay fraction [%w]
snippet = "OpenLandMap/SOL/SOL_CLAY-WFRACTION_USDA-3A1A1A_M/v02"
# Define the scale factor in accordance with the dataset description.
scale_factor = 1 * 0.01

elif param == "orgc":  # Organic Carbon fraction [g/kg]
snippet = "OpenLandMap/SOL/SOL_ORGANIC-CARBON_USDA-6A1C_M/v02"
# Define the scale factor in accordance with the dataset description.
scale_factor = 5 * 0.001  # to get kg/kg
else:
return print("error")

# Apply the scale factor to the ee.Image.
dataset = ee.Image(snippet).multiply(scale_factor)

return dataset

We apply this function to import soil properties:

# Image associated with the sand content.
sand = get_soil_prop("sand")

# Image associated with the clay content.
clay = get_soil_prop("clay")

# Image associated with the organic carbon content.
orgc = get_soil_prop("orgc")

To illustrate the result, we define a new method for handing Earth Engine tiles and using it to display the clay content of the soil at a given reference depth, to a Leaflet map.

def add_ee_layer(self, ee_image_object, vis_params, name):
"""Adds a method for displaying Earth Engine image tiles to folium map."""
map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)
folium.raster_layers.TileLayer(
tiles=map_id_dict["tile_fetcher"].url_format,
name=name,
overlay=True,
control=True,

# Add Earth Engine drawing method to folium.

my_map = folium.Map(location=[lat, lon], zoom_start=3)

# Set visualization parameters.
vis_params = {
"bands": ["b0"],
"min": 0.01,
"max": 1,
"opacity": 1,
"palette": ["white", "#464646"],
}

# Add the sand content data to the map object.

# Add a marker at the location of interest.

# Add a layer control panel to the map.

# Display the map.
display(my_map)

Now, a function is defined to get soil properties at a given location. The following function returns a dictionary indicating the value of the parameter of interest for each standard depth (in centimeter). This function uses the ee.Image.sample method to evaluate the ee.Image properties on the region of interest. The result is then transferred client-side using the ee.Image.getInfo method.

In the example below, we are asking for the sand content.

def local_profile(dataset, poi, buffer):
# Get properties at the location of interest and transfer to client-side.
prop = dataset.sample(poi, buffer).select(olm_bands).getInfo()

# Selection of the features/properties of interest.
profile = prop["features"][0]["properties"]

# Re-shaping of the dict.
profile = {key: round(val, 3) for key, val in profile.items()}

return profile

# Apply the function to get the sand profile.
profile_sand = local_profile(sand, poi, scale)

# Print the result.
print("Sand content profile at the location of interest:\n", profile_sand)

Sand content profile at the location of interest:
{'b0': 0.36, 'b10': 0.35, 'b100': 0.37, 'b200': 0.39, 'b30': 0.34, 'b60': 0.35}

We now apply the function to plot the profile of the soil regarding sand and clay and organic carbon content at the location of interest:

# Clay and organic content profiles.
profile_clay = local_profile(clay, poi, scale)
profile_orgc = local_profile(orgc, poi, scale)

# Data visualization in the form of a bar plot.
fig, ax = plt.subplots(figsize=(15, 6))
ax.axes.get_yaxis().set_visible(False)

# Definition of label locations.
x = np.arange(len(olm_bands))

# Definition of the bar width.
width = 0.25

# Bar plot representing the sand content profile.
rect1 = ax.bar(
x - width,
[round(100 * profile_sand[b], 2) for b in olm_bands],
width,
label="Sand",
color="#ecebbd",
)

# Bar plot representing the clay content profile.
rect2 = ax.bar(
x,
[round(100 * profile_clay[b], 2) for b in olm_bands],
width,
label="Clay",
color="#6f6c5d",
)

# Bar plot representing the organic carbon content profile.
rect3 = ax.bar(
x + width,
[round(100 * profile_orgc[b], 2) for b in olm_bands],
width,
label="Organic Carbon",
color="black",
alpha=0.75,
)

# Definition of a function to attach a label to each bar.
def autolabel_soil_prop(rects):
"""Attach a text label above each bar in *rects*, displaying its height."""
for rect in rects:
height = rect.get_height()
ax.annotate(
"{}".format(height) + "%",
xy=(rect.get_x() + rect.get_width() / 2, height),
xytext=(0, 3),  # 3 points vertical offset.
textcoords="offset points",
ha="center",
va="bottom",
fontsize=10,
)

# Application of the function to each barplot.
autolabel_soil_prop(rect1)
autolabel_soil_prop(rect2)
autolabel_soil_prop(rect3)

# Title of the plot.
ax.set_title("Properties of the soil at different depths (mass content)", fontsize=14)

# Properties of x/y labels and ticks.
ax.set_xticks(x)
x_labels = [str(d) + " cm" for d in olm_depths]
ax.set_xticklabels(x_labels, rotation=45, fontsize=10)

ax.spines["left"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)

# Shrink current axis's height by 10% on the bottom.
box = ax.get_position()
ax.set_position([box.x0, box.y0 + box.height * 0.1, box.width, box.height * 0.9])

# Add a legend below current axis.
ax.legend(
loc="upper center", bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=3
)

plt.show()

#### Expression to calculate hydraulic properties

Now that soil properties are described, the water content at the field capacity and at the wilting point can be calculated according to the equation defined at the beginning of this section. Please note that in the equation of Saxton & Rawls (2006), the wilting point and field capacity are calculated using the Organic Matter content ($$OM$$) and not the Organic Carbon content ($$OC$$). In the following, we convert $$OC$$ into $$OM$$ using the corrective factor known as the Van Bemmelen factor:

$$0M = 1.724 \times OC$$

Several operators are available to perform basic mathematical operations on image bands: add(), subtract(), multiply() and divide(). Here, we multiply the organic content by the Van Bemmelen factor. It is done using the ee.Image.multiply method on the organic carbon content ee.Image.

# Conversion of organic carbon content into organic matter content.
orgm = orgc.multiply(1.724)

# Organic matter content profile.
profile_orgm = local_profile(orgm, poi, scale)
print("Organic Matter content profile at the location of interest:\n", profile_orgm)

Organic Matter content profile at the location of interest:
{'b0': 0.043, 'b10': 0.043, 'b100': 0.009, 'b200': 0.009, 'b30': 0.017, 'b60': 0.009}

When the mathematical operation to apply to the ee.Image becomes too complex, the ee.Image.expression is a good alternative. We use it in the following code block since the calculation of wilting point and field capacity relies on multiple parameters and images. This method takes two arguments:

• a string formalizing the arithmetic expression we want to evaluate,
• a dict associating images to each parameter of the arithmetic expression.

The mathematical expression is applied as follows to determine wilting point and field capacity:

# Initialization of two constant images for wilting point and field capacity.
wilting_point = ee.Image(0)
field_capacity = ee.Image(0)

# Calculation for each standard depth using a loop.
for key in olm_bands:
# Getting sand, clay and organic matter at the appropriate depth.
si = sand.select(key)
ci = clay.select(key)
oi = orgm.select(key)

# Calculation of the wilting point.
# The theta_1500t parameter is needed for the given depth.
theta_1500ti = (
ee.Image(0)
.expression(
"-0.024 * S + 0.487 * C + 0.006 * OM + 0.005 * (S * OM)\
- 0.013 * (C * OM) + 0.068 * (S * C) + 0.031",
{
"S": si,
"C": ci,
"OM": oi,
},
)
.rename("T1500ti")
)

# Final expression for the wilting point.
wpi = theta_1500ti.expression(
"T1500ti + ( 0.14 * T1500ti - 0.002)", {"T1500ti": theta_1500ti}
).rename("wpi")

# Add as a new band of the global wilting point ee.Image.
# Do not forget to cast the type with float().

# Same process for the calculation of the field capacity.
# The parameter theta_33t is needed for the given depth.
theta_33ti = (
ee.Image(0)
.expression(
"-0.251 * S + 0.195 * C + 0.011 * OM +\
0.006 * (S * OM) - 0.027 * (C * OM)+\
0.452 * (S * C) + 0.299",
{
"S": si,
"C": ci,
"OM": oi,
},
)
.rename("T33ti")
)

# Final expression for the field capacity of the soil.
fci = theta_33ti.expression(
"T33ti + (1.283 * T33ti * T33ti - 0.374 * T33ti - 0.015)",
{"T33ti": theta_33ti.select("T33ti")},
)

# Add a new band of the global field capacity ee.Image.

Let's see the result around our location of interest:

profile_wp = local_profile(wilting_point, poi, scale)
profile_fc = local_profile(field_capacity, poi, scale)

print("Wilting point profile:\n", profile_wp)
print("Field capacity profile:\n", profile_fc)

Wilting point profile:
{'b0': 0.152, 'b10': 0.158, 'b100': 0.175, 'b200': 0.169, 'b30': 0.175, 'b60': 0.181}
Field capacity profile:
{'b0': 0.271, 'b10': 0.278, 'b100': 0.289, 'b200': 0.28, 'b30': 0.294, 'b60': 0.297}

The result is displayed using barplots as follows:

fig, ax = plt.subplots(figsize=(15, 6))
ax.axes.get_yaxis().set_visible(False)

# Definition of the label locations.
x = np.arange(len(olm_bands))

# Width of the bar of the barplot.
width = 0.25

# Barplot associated with the water content at the wilting point.
rect1 = ax.bar(
x - width / 2,
[round(profile_wp[b] * 100, 2) for b in olm_bands],
width,
label="Water content at wilting point",
color="red",
alpha=0.5,
)

# Barplot associated with the water content at the field capacity.
rect2 = ax.bar(
x + width / 2,
[round(profile_fc[b] * 100, 2) for b in olm_bands],
width,
label="Water content at field capacity",
color="blue",
alpha=0.5,
)

# Add Labels on top of bars.
autolabel_soil_prop(rect1)
autolabel_soil_prop(rect2)

# Title of the plot.
ax.set_title("Hydraulic properties of the soil at different depths", fontsize=14)

# Properties of x/y labels and ticks.
ax.set_xticks(x)
x_labels = [str(d) + " cm" for d in olm_depths]
ax.set_xticklabels(x_labels, rotation=45, fontsize=10)

ax.spines["left"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)

# Shrink current axis's height by 10% on the bottom.
box = ax.get_position()
ax.set_position([box.x0, box.y0 + box.height * 0.1, box.width, box.height * 0.9])

# Put a legend below current axis.
ax.legend(
loc="upper center", bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=2
)

plt.show()

### Getting meteorological datasets

#### Datasets exploration

The meteorological data used in our implementation of the TM procedure relies on the following datasets:

Both datasets are imported as follows, specifying the bands of interest using .select() and the period of interest using .filterDate().

# Import precipitation.
pr = (
ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
.select("precipitation")
.filterDate(i_date, f_date)
)

# Import potential evaporation PET and its quality indicator ET_QC.
pet = (
ee.ImageCollection("MODIS/006/MOD16A2")
.select(["PET", "ET_QC"])
.filterDate(i_date, f_date)
)

Now we can have a closer look around our location of interest. To evaluate the properties of an ee.ImageCollection, the ee.ImageCollection.getRegion method is used and combined with ee.ImageCollection.getInfo method for a client-side visualization.

# Evaluate local precipitation conditions.
local_pr = pr.getRegion(poi, scale).getInfo()
pprint.pprint(local_pr[:5])

[['id', 'longitude', 'latitude', 'time', 'precipitation'],
['20150101', 5.142855001584261, 45.77365530231022, 1420070400000, 0],
['20150102', 5.142855001584261, 45.77365530231022, 1420156800000, 0],
['20150103', 5.142855001584261, 45.77365530231022, 1420243200000, 0],
['20150104', 5.142855001584261, 45.77365530231022, 1420329600000, 0]]

We now establish a procedure to get meteorological data around a given location in form of a pandas.DataFrame:

def ee_array_to_df(arr, list_of_bands):
"""Transforms client-side ee.Image.getRegion array to pandas.DataFrame."""
df = pd.DataFrame(arr)

# Convert the data to numeric values.
for band in list_of_bands:
df[band] = pd.to_numeric(df[band], errors="coerce")

# Convert the time field into a datetime.
df["datetime"] = pd.to_datetime(df["time"], unit="ms")

# Keep the columns of interest.
df = df[["time", "datetime", *list_of_bands]]

# The datetime column is defined as index.
df = df.set_index("datetime")

return df

We apply the function and see the head of the resulting pandas.DataFrame:

pr_df = ee_array_to_df(local_pr, ["precipitation"])

We do the same for potential evaporation:

# Evaluate local potential evapotranspiration.
local_pet = pet.getRegion(poi, scale).getInfo()

# Transform the result into a pandas dataframe.
pet_df = ee_array_to_df(local_pet, ["PET", "ET_QC"])

Looking at both pandas.DataFrame shows the following points:

• the time resolution between both datasets is not the same,
• for some reasons, potential evapotranspiration cannot be calculated at some dates. It corresponds to lines where the quality indicator ET_QC is higher than 1.

Both issues must be handled before implementing the iterative process: we want to work on a similar timeline with potential evapotranspiration and precipitation, and we want to avoid missing values.

#### Resampling the time resolution of an ee.ImageCollection

To address these issues (homogeneous time index and missing values), we make a sum resampling of both datasets by month. When PET cannot be calculated, the monthly averaged value is considered. The key steps and functions used to resample are described below:

• A new date index is defined as a sequence using the ee.List.sequence method.
• A function representing the resampling operation is defined. This function consists of grouping images of the desired time interval and calculating the sum. The sum is calculated by taking the mean between available images and multiplying it by the duration of the interval.
• The user-supplied function is then mapped over the new time index using .map().

Finally, the resampling procedure reads as follows:

def sum_resampler(coll, freq, unit, scale_factor, band_name):
"""
This function aims to resample the time scale of an ee.ImageCollection.
The function returns an ee.ImageCollection with the averaged sum of the
band on the selected frequency.

coll: (ee.ImageCollection) only one band can be handled
freq: (int) corresponds to the resampling frequence
unit: (str) corresponds to the resampling time unit.
must be 'day', 'month' or 'year'
scale_factor (float): scaling factor used to get our value in the good unit
band_name (str) name of the output band
"""

# Define initial and final dates of the collection.
firstdate = ee.Date(
coll.sort("system:time_start", True).first().get("system:time_start")
)

lastdate = ee.Date(
coll.sort("system:time_start", False).first().get("system:time_start")
)

# Calculate the time difference between both dates.
diff_dates = lastdate.difference(firstdate, unit)

# Define a new time index (for output).
new_index = ee.List.sequence(0, ee.Number(diff_dates), freq)

# Define the function that will be applied to our new time index.
def apply_resampling(date_index):
# Define the starting date to take into account.

# Define the ending date to take into account according
# to the desired frequency.

# Calculate the number of days between starting and ending days.
diff_days = enddate.difference(startdate, "day")

# Calculate the composite image.
image = (
coll.filterDate(startdate, enddate)
.mean()
.multiply(diff_days)
.multiply(scale_factor)
.rename(band_name)
)

# Return the final image with the appropriate time index.
return image.set("system:time_start", startdate.millis())

# Map the function to the new time index.
res = new_index.map(apply_resampling)

# Transform the result into an ee.ImageCollection.
res = ee.ImageCollection(res)

return res

The precipitation dataset is now resampled by month as follows:

• the collection to resample is defined as pr,
• we want a collection on a monthly basis then freq = 1 and unit = "month",
• there is no correction factor to apply according to the dataset description then scale_factor = 1,
• "pr" is the name of the output band.
# Apply the resampling function to the precipitation dataset.
pr_m = sum_resampler(pr, 1, "month", 1, "pr")

# Evaluate the result at the location of interest.
pprint.pprint(pr_m.getRegion(poi, scale).getInfo()[:5])

[['id', 'longitude', 'latitude', 'time', 'pr'],
['0', 5.142855001584261, 45.77365530231022, 1420070400000, 83.12794733047485],
['1', 5.142855001584261, 45.77365530231022, 1422748800000, 61.235233306884766],
['2', 5.142855001584261, 45.77365530231022, 1425168000000, 60.416391015052795],
['3', 5.142855001584261, 45.77365530231022, 1427846400000, 48.404173851013184]]

For evapotranspiration, we have to be careful with the unit. The dataset gives us an 8-day sum and a scale factor of 10 is applied. Then, to get a homogeneous unit, we need to rescale by dividing by 8 and 10: $$\frac{1}{10 \times 8} = 0.0125$$.

# Apply the resampling function to the PET dataset.
pet_m = sum_resampler(pet.select("PET"), 1, "month", 0.0125, "pet")

# Evaluate the result at the location of interest.
pprint.pprint(pet_m.getRegion(poi, scale).getInfo()[:5])

[['id', 'longitude', 'latitude', 'time', 'pet'],
['0', 5.142855001584261, 45.77365530231022, 1420070400000, 23.637500000000003],
['1', 5.142855001584261, 45.77365530231022, 1422748800000, 39.2],
['2', 5.142855001584261, 45.77365530231022, 1425168000000, 76.434375],
['3', 5.142855001584261, 45.77365530231022, 1427846400000, 137]]

We now combine both ee.ImageCollection objects (pet_m and pr_m) using the ee.ImageCollection.combine method. Note that corresponding images in both ee.ImageCollection objects need to have the same time index before combining.

# Combine precipitation and evapotranspiration.
meteo = pr_m.combine(pet_m)

# Import meteorological data as an array at the location of interest.
meteo_arr = meteo.getRegion(poi, scale).getInfo()

# Print the result.
pprint.pprint(meteo_arr[:5])

[['id', 'longitude', 'latitude', 'time', 'pr', 'pet'],
['0',
5.142855001584261,
45.77365530231022,
1420070400000,
83.12794733047485,
23.637500000000003],
['1',
5.142855001584261,
45.77365530231022,
1422748800000,
61.235233306884766,
39.2],
['10',
5.142855001584261,
45.77365530231022,
1446336000000,
35.43960928916931,
39.5],
['11',
5.142855001584261,
45.77365530231022,
1448928000000,
16.603714764118195,
25.058333333333337]]

We evaluate the result on our location of interest:

# Transform the array into a pandas dataframe and sort the index.
meteo_df = ee_array_to_df(meteo_arr, ["pr", "pet"]).sort_index()

# Data visualization
fig, ax = plt.subplots(figsize=(15, 6))

# Barplot associated with precipitations.
meteo_df["pr"].plot(kind="bar", ax=ax, label="precipitation")

# Barplot associated with potential evapotranspiration.
meteo_df["pet"].plot(
kind="bar", ax=ax, label="potential evapotranspiration", color="orange", alpha=0.5
)

ax.legend()

ax.set_ylabel("Intensity [mm]")
ax.set_xlabel(None)

# Define the date format and shape of x-labels.
x_labels = meteo_df.index.strftime("%m-%Y")
ax.set_xticklabels(x_labels, rotation=90, fontsize=10)

plt.show()

### Implementation of the TM procedure

#### Description

Some additional definitions are needed to formalize the Thornthwaite-Mather procedure. The following definitions are given in accordance with Allen et al. (1998) (the document can be downloaded here):

$$TAW = 1000 \times (\theta_{FC} − \theta_{WP})\times Z{r}$$

where:

• $$TAW$$: the total available soil water in the root zone [$$mm$$],
• $$\theta_{FC}$$: the water content at the field capacity [$$m^{3} m^{-3}$$],
• $$\theta_{WP}$$: the water content at wilting point [$$m^{3} m^{-3}$$],
• $$Z_{r}$$: the rooting depth [$$mm$$],

Typical values of $$\theta_{FC}$$ and $$\theta_{WP}$$ for different soil types are given in Table 19 of Allen et al. (1998).

The readily available water ($$RAW$$) is given by $$RAW = p \times TAW$$, where $$p$$ is the average fraction of $$TAW$$ that can be depleted from the root zone before moisture stress (ranging between 0 to 1). This quantity is also noted $$ST_{FC}$$ which is the available water stored at field capacity in the root zone.

Ranges of maximum effective rooting depth $$Z_{r}$$, and soil water depletion fraction for no stress $$p$$, for common crops are given in the Table 22 of Allen et al. (1998). In addition, a global effective plant rooting depth dataset is provided by Yang et al. (2016) with a resolution of 0.5° by 0.5° (see the paper here and the dataset here).

According to this global dataset, the effective rooting depth around our region of interest (France) can reasonably assumed to $$Z_{r} = 0.5$$. Additionally, the parameter $$p$$ is also assumed constant and equal to and $$p = 0.5$$ which is in line with common values described in Table 22 of Allen et al. (1998).

zr = ee.Image(0.5)
p = ee.Image(0.5)

In the following, we also consider an averaged value between reference depths of the water content at wilting point and field capacity:

def olm_prop_mean(olm_image, band_output_name):
"""
This function calculates an averaged value of
soil properties between reference depths.
"""
mean_image = olm_image.expression(
"(b0 + b10 + b30 + b60 + b100 + b200) / 6",
{
"b0": olm_image.select("b0"),
"b10": olm_image.select("b10"),
"b30": olm_image.select("b30"),
"b60": olm_image.select("b60"),
"b100": olm_image.select("b100"),
"b200": olm_image.select("b200"),
},
).rename(band_output_name)

return mean_image

# Apply the function to field capacity and wilting point.
fcm = olm_prop_mean(field_capacity, "fc_mean")
wpm = olm_prop_mean(wilting_point, "wp_mean")

# Calculate the theoretical available water.
taw = (
(fcm.select("fc_mean").subtract(wpm.select("wp_mean"))).multiply(1000).multiply(zr)
)

# Calculate the stored water at the field capacity.
stfc = taw.multiply(p)

The Thornthwaite-Mather procedure used to estimate groundwater recharge is explicitly described by Steenhuis and Van der Molen (1985). This procedure uses monthly sums of potential evaporation, cumulative precipitation, and the moisture status of the soil which is calculated iteratively. The moisture status of the soils depends on the accumulated potential water loss ($$APWL$$). This parameter is calculated depending on whether the potential evaporation is greater than or less than the cumulative precipitation. The procedure reads as follow:

Case 1: potential evapotranspiration is higher than precipitation.

In that case, $$PET>P$$ and $$APWL_{m}$$ is incremented as follows: $$APWL_{m} = APWL_{m - 1} + (PET_{m} - P_{m})$$ where:

• $$APWL_{m}$$ (respectively $$APWL_{m - 1}$$) represents the accumulated potential water loss for the month $$m$$ (respectively at the previous month $$m - 1$$)
• $$PET_{m}$$ the cumulative potential evapotranspiration at month $$m$$,
• $$P_{m}$$ the cumulative precipitation at month $$m$$,

and the relationship between $$APWL$$ and the amount of water stored in the root zone for the month $$m$$ is expressed as: $$ST_{m} = ST_{FC} \times [\textrm{exp}(-APWL_{m}/ST_{FC})]$$ where $$ST_{m}$$ is the available water stored in the root zone for the month $$m$$.

Case 2: potential evapotranspiration is lower than precipitation.

In that case, $$PET<P$$ and $$ST_{m}$$ is incremented as follows: $$ST_{m} = ST_{m-1} + (P_{m} - PET_{m})$$.

Case 2.1: the storage $$ST_{m}$$ is higher than the water stored at the field capacity.

If $$ST_{m} > ST_{FC}$$ the recharge is calculated as: $$R_{m} = ST_{m} - ST_{FC} + P_{m} - PET_{m}$$

In addition, the water stored at the end of the month $$m$$ becomes equal to $$ST_{FC}$$ and $$APWL_{m}$$ is set equal to zero.

Case 2.2: the storage $$ST_{m}$$ is less than or equal to the water stored at the field capacity.

If $$ST_{m} <= ST_{FC}$$, $$APWL_{m}$$ is implemented as follows: $$APWL_{m} = ST_{FC} \times \textrm{ln}(ST_{m}/ST_{FC})$$, and no percolation occurs.

#### Initialization

The initial time of the calculation is defined according to the first date of the meteorological collection:

# Define the initial time (time0) according to the start of the collection.
time0 = meteo.first().get("system:time_start")

Then, we initialize the calculation with an ee.Image where all bands associated to the hydric state of the soil are set equal to ee.Image(0), except for the initial storage which is considered to be equal to the water content at field capacity, meaning that $$ST_{0} = ST_{FC}$$.

# Initialize all bands describing the hydric state of the soil.
# Do not forget to cast the type of the data with a .float().
# Initial recharge.
initial_rech = ee.Image(0).set("system:time_start", time0).select([0], ["rech"]).float()

# Initialization of APWL.
initial_apwl = ee.Image(0).set("system:time_start", time0).select([0], ["apwl"]).float()

# Initialization of ST.
initial_st = stfc.set("system:time_start", time0).select([0], ["st"]).float()

# Initialization of precipitation.
initial_pr = ee.Image(0).set("system:time_start", time0).select([0], ["pr"]).float()

# Initialization of potential evapotranspiration.
initial_pet = ee.Image(0).set("system:time_start", time0).select([0], ["pet"]).float()

We combine all these bands into one ee.Image adding new bands to the first using the ee.Image.addBands method:

initial_image = initial_rech.addBands(
ee.Image([initial_apwl, initial_st, initial_pr, initial_pet])
)

We also initialize a list in which new images will be added after each iteration. We create this server-side list using the ee.List method.

image_list = ee.List([initial_image])

#### Iteration over an ee.ImageCollection

The procedure is implemented by means of the ee.ImageCollection.iterate method, which applies a user-supplied function to each element of a collection. For each time step, groundwater recharge is calculated using the recharge_calculator considering the previous hydric state of the soil and current meteorological conditions.

Of course, considering the TM description, several cases must be distinguished to calculate groundwater recharge. The distinction is made by the definition of binary layers with different logical operations. It allows specific calculations to be applied in areas where a given condition is true using the ee.Image.where method.

The function we apply to each element of the meteorological dataset to calculate groundwater recharge is defined as follows.

def recharge_calculator(image, image_list):
"""
Contains operations made at each iteration.
"""
# Determine the date of the current ee.Image of the collection.
localdate = image.date().millis()

# Import previous image stored in the list.
prev_im = ee.Image(ee.List(image_list).get(-1))

# Import previous APWL and ST.
prev_apwl = prev_im.select("apwl")
prev_st = prev_im.select("st")

# Import current precipitation and evapotranspiration.
pr_im = image.select("pr")
pet_im = image.select("pet")

# Initialize the new bands associated with recharge, apwl and st.
# DO NOT FORGET TO CAST THE TYPE WITH .float().
new_rech = (
ee.Image(0)
.set("system:time_start", localdate)
.select([0], ["rech"])
.float()
)

new_apwl = (
ee.Image(0)
.set("system:time_start", localdate)
.select([0], ["apwl"])
.float()
)

new_st = (
prev_st.set("system:time_start", localdate).select([0], ["st"]).float()
)

# Calculate bands depending on the situation using binary layers with
# logical operations.

# CASE 1.
# Define zone1: the area where PET > P.
zone1 = pet_im.gt(pr_im)

# Calculation of APWL in zone 1.
# Implementation of zone 1 values for APWL.
new_apwl = new_apwl.where(zone1, zone1_apwl)

# Calculate ST in zone 1.
zone1_st = prev_st.multiply(
ee.Image.exp(zone1_apwl.divide(stfc).multiply(-1))
).rename("st")
# Implement ST in zone 1.
new_st = new_st.where(zone1, zone1_st)

# CASE 2.
# Define zone2: the area where PET <= P.
zone2 = pet_im.lte(pr_im)

# Calculate ST in zone 2.
# Implement ST in zone 2.
new_st = new_st.where(zone2, zone2_st)

# CASE 2.1.
# Define zone21: the area where PET <= P and ST >= STfc.
zone21 = zone2.And(zone2_st.gte(stfc))

# Calculate recharge in zone 21.
zone21_re = zone2_st.subtract(stfc).rename("rech")
# Implement recharge in zone 21.
new_rech = new_rech.where(zone21, zone21_re)
# Implement ST in zone 21.
new_st = new_st.where(zone21, stfc)

# CASE 2.2.
# Define zone 22: the area where PET <= P and ST < STfc.
zone22 = zone2.And(zone2_st.lt(stfc))

# Calculate APWL in zone 22.
zone22_apwl = (
stfc.multiply(-1).multiply(ee.Image.log(zone2_st.divide(stfc))).rename("apwl")
)
# Implement APWL in zone 22.
new_apwl = new_apwl.where(zone22, zone22_apwl)

# Create a mask around area where recharge can effectively be calculated.
# Where we have have PET, P, FCm, WPm (except urban areas, etc.).

# Add all Bands to our ee.Image.
new_image = new_rech.addBands(ee.Image([new_apwl, new_st, pr_im, pet_im]))

# Add the new ee.Image to the ee.List.

The TM procedure can now be applied to the meteorological ee.ImageCollection:

# Iterate the user-supplied function to the meteo collection.
rech_list = meteo.iterate(recharge_calculator, image_list)

# Remove the initial image from our list.
rech_list = ee.List(rech_list).remove(initial_image)

# Transform the list into an ee.ImageCollection.
rech_coll = ee.ImageCollection(rech_list)

Let's have a look at the result around the location of interest:

arr = rech_coll.getRegion(poi, scale).getInfo()
rdf = ee_array_to_df(arr, ["pr", "pet", "apwl", "st", "rech"]).sort_index()

The result can be displayed in the form of a barplot as follows:

# Data visualization in the form of barplots.
fig, ax = plt.subplots(figsize=(15, 6))

# Barplot associated with precipitation.
rdf["pr"].plot(kind="bar", ax=ax, label="precipitation", alpha=0.5)

# Barplot associated with potential evapotranspiration.
rdf["pet"].plot(
kind="bar", ax=ax, label="potential evapotranspiration", color="orange", alpha=0.2
)

# Barplot associated with groundwater recharge
rdf["rech"].plot(kind="bar", ax=ax, label="recharge", color="green", alpha=1)

ax.legend()

# Define x/y-labels properties.
ax.set_ylabel("Intensity [mm]")
ax.set_xlabel(None)

# Define the date format and shape of x-labels.
x_labels = rdf.index.strftime("%m-%Y")
ax.set_xticklabels(x_labels, rotation=90, fontsize=10)

plt.show()

The result shows the distribution of precipitation, potential evapotranspiration, and groundwater recharge along the year. It shows that in our area of interest, groundwater recharge generally occurs from October to March. Even though a significant amount of precipitation occurs from April to September, evapotranspiration largely dominates because of high temperatures and sun exposure during these months. The result is that no percolation into aquifers occurs during this period.

Now the annual average recharge over the period of interest can be calculated. To do that, we resample the DataFrame we've just created:

# Resample the pandas dataframe on a yearly basis making the sum by year.
rdfy = rdf.resample("Y").sum()

# Calculate the mean value.
mean_recharge = rdfy["rech"].mean()

# Print the result.
print(
"The mean annual recharge at our point of interest is", int(mean_recharge), "mm/an"
)

The mean annual recharge at our point of interest is 154 mm/an

#### Groundwater recharge comparison between multiple places

We now may want to get local information about groundwater recharge and/or map this variable on an area of interest.

Let's define a function to get the local information based on the ee.ImageCollection we've just built:

def get_local_recharge(i_date, f_date, lon, lat, scale):
"""
Returns a pandas df describing the cumulative groundwater
recharge by month
"""
# Define the location of interest with a point.
poi = ee.Geometry.Point(lon, lat)

# Evaluate the recharge around the location of interest.
rarr = rech_coll.filterDate(i_date, f_date).getRegion(poi, scale).getInfo()

# Transform the result into a pandas dataframe.
rdf = ee_array_to_df(rarr, ["pr", "pet", "apwl", "st", "rech"]).sort_index()
return rdf

We now use this function on a second point of interest located near the city of Montpellier (France). This city is located in the south of France, and precipitation and groundwater recharge are expected to be much lower than in the previous case.

# Define the second location of interest by longitude/latitude.
lon2 = 4.137152
lat2 = 43.626945

# Calculate the local recharge condition at this location.
rdf2 = get_local_recharge(i_date, f_date, lon2, lat2, scale)

# Resample the resulting pandas dataframe on a yearly basis (sum by year).
rdf2y = rdf2.resample("Y").sum()

# Data Visualization
fig, ax = plt.subplots(figsize=(15, 6))
ax.axes.get_yaxis().set_visible(False)

# Define the x-label locations.
x = np.arange(len(rdfy))

# Define the bar width.
width = 0.25

# Bar plot associated to groundwater recharge at the 1st location of interest.
rect1 = ax.bar(
x - width / 2, rdfy.rech, width, label="Lyon (France)", color="blue", alpha=0.5
)

# Bar plot associated to groundwater recharge at the 2nd location of interest.
rect2 = ax.bar(
x + width / 2,
rdf2y.rech,
width,
label="Montpellier (France)",
color="red",
alpha=0.5,
)

# Define a function to attach a label to each bar.
def autolabel_recharge(rects):
"""Attach a text label above each bar in *rects*, displaying its height."""
for rect in rects:
height = rect.get_height()
ax.annotate(
"{}".format(int(height)) + " mm",
xy=(rect.get_x() + rect.get_width() / 2, height),
xytext=(0, 3),  # 3 points vertical offset
textcoords="offset points",
ha="center",
va="bottom",
fontsize=8,
)

autolabel_recharge(rect1)
autolabel_recharge(rect2)

# Calculate the averaged annual recharge at both locations of interest.
place1mean = int(rdfy["rech"].mean())
place2mean = int(rdf2y["rech"].mean())

# Add an horizontal line associated with averaged annual values (location 1).
ax.hlines(
place1mean,
xmin=min(x) - width,
xmax=max(x) + width,
color="blue",
lw=0.5,
label="average " + str(place1mean) + " mm/y",
alpha=0.5,
)

# Add an horizontal line associated with averaged annual values (location 2).
ax.hlines(
place2mean,
xmin=min(x) - width,
xmax=max(x) + width,
color="red",
lw=0.5,
label="average " + str(place2mean) + " mm/y",
alpha=0.5,
)

ax.set_title("Groundwater recharge comparison between two places", fontsize=12)

# Define some x/y-axis properties.
ax.set_xticks(x)
x_labels = rdfy.index.year.tolist()
ax.set_xticklabels(x_labels, rotation=45, fontsize=10)

ax.spines["left"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)

# Shrink current axis's height by 10% on the bottom.
box = ax.get_position()
ax.set_position([box.x0, box.y0 + box.height * 0.1, box.width, box.height * 0.9])

# Add a legend below current axis.
ax.legend(
loc="upper center", bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=2
)

plt.show()

The result shows that the annual recharge in Lyon is almost twice as high as in the area of Montpellier. The result also shows a great variability of the annual recharge ranging from 98 mm/y to 258 mm/y in Lyon and from 16 mm/y to 147 mm/y in Montpellier.

### Groundwater recharge map of France

To get a map of groundwater recharge around our region of interest, let's create a mean composite ee.Image based on our resulting ee.ImageCollection.

# Calculate the averaged annual recharge.
annual_rech = rech_coll.select("rech").mean().multiply(12)

# Calculate the average annual precipitation.
annual_pr = rech_coll.select("pr").mean().multiply(12)

# Get a feature collection of administrative boundaries.

# Filter the feature collection to subset France.

# Clip the composite ee.Images around the region of interest.
rech_france = annual_rech.clip(france)
pr_france = annual_pr.clip(france)

And finally the map can be drawn.

# Create a folium map.
my_map = folium.Map(location=[lat, lon], zoom_start=6, zoom_control=False)

# Set visualization parameters for recharge.
rech_vis_params = {
"bands": "rech",
"min": 0,
"max": 500,
"opacity": 1,
"palette": ["red", "orange", "yellow", "green", "blue", "purple"],
}

# Set visualization parameters for precipitation.
pr_vis_params = {
"bands": "pr",
"min": 500,
"max": 1500,
"opacity": 1,
"palette": ["white", "blue"],
}

# Define a recharge colormap.
rech_colormap = cm.LinearColormap(
colors=rech_vis_params["palette"],
vmin=rech_vis_params["min"],
vmax=rech_vis_params["max"],
)

# Define a precipitation colormap.
pr_colormap = cm.LinearColormap(
colors=pr_vis_params["palette"],
vmin=pr_vis_params["min"],
vmax=pr_vis_params["max"],
)

# Caption of the recharge colormap.
rech_colormap.caption = "Average annual recharge rate (mm/year)"

# Caption of the precipitation colormap.
pr_colormap.caption = "Average annual precipitation rate (mm/year)"

# Add the precipitation composite to the map object.

# Add the recharge composite to the map object.

# Add a marker at both locations of interest.

# Add the colormaps to the map.

# Add a layer control panel to the map.

# Display the map.
display(my_map)

The resulting map represents the averaged annual precipitation and groundwater recharge rate in France over the period 2015-2019. Groundwater recharge is represented between 0 (red) to 500 mm/year (purple). Precipitation is represented between 500 (white) to 1500 mm/year (blue).

• The precipitation dataset used in this tutorial is not available over +/- 50 deg latitudes. That is why (1) the northern part of France is not covered by the precipitation/recharge, and (2) some anomalies/discontinuities can be observed near the boundary of the dataset.
• The procedure described in this tutorial does not consider the terrain slope which has a strong influence on runoff, especially in areas with significant relief. In these areas, groundwater recharge is overestimated because runoff is not considered.
• Snow cover has a more complex contribution to groundwater recharge than precipitation. This process is not taken into account in this tutorial. Hence, groundwater recharge is overestimated in areas where snow cover period is significant.
• The recharge calculation can be highly influenced by missing precipitation or potential evapotranspiration values. This issue is partially addressed by the dataset resampling procedure, but missing value can still exist and this point must be checked before considering that the calculated recharge is representative and reliable.

### Bibliography

Allen RG, Pereira LS, Raes D, Smith M (1998). Crop evapotranspiration: guidelines for computing crop water requirements. Irrigation and Drainage Paper 56, FAO, Rome.

Saxton, K. E., & Rawls, W. J. (2006). Soil water characteristic estimates by texture and organic matter for hydrologic solutions. Soil science society of America Journal, 70(5), 1569-1578.

Steenhuis, T. S., & Van der Molen, W. H. (1986). The Thornthwaite-Mather procedure as a simple engineering method to predict recharge. Journal of Hydrology, 84(3-4), 221-229.

Thornthwaite, C. W., & Mather, J. R. (1957). Instructions and tables for computing potential evapotranspiration and the water balance. Publ. Climatol., 10(3).

Yang, Y., Donohue, R. J., & McVicar, T. R. (2016). Global estimation of effective plant rooting depth: Implications for hydrological modeling. Water Resources Research, 52(10), 8260-8276.

### Acknowledgements

Thanks to Susanne Benz for precious advice.

[{ "type": "thumb-down", "id": "missingTheInformationINeed", "label":"Missing the information I need" },{ "type": "thumb-down", "id": "tooComplicatedTooManySteps", "label":"Too complicated / too many steps" },{ "type": "thumb-down", "id": "outOfDate", "label":"Out of date" },{ "type": "thumb-down", "id": "samplesCodeIssue", "label":"Samples / code issue" },{ "type": "thumb-down", "id": "otherDown", "label":"Other" }]
[{ "type": "thumb-up", "id": "easyToUnderstand", "label":"Easy to understand" },{ "type": "thumb-up", "id": "solvedMyProblem", "label":"Solved my problem" },{ "type": "thumb-up", "id": "otherUp", "label":"Other" }]


text:
# §28.11 Expansions in Series of Mathieu Functions

Let $f(z)$ be a $2\pi$-periodic function that is analytic in an open doubly-infinite strip $S$ that contains the real axis, and $q$ be a normal value (§28.7). Then

 28.11.1 $f(z)=\alpha_{0}\mathop{\mathrm{ce}_{0}\/}\nolimits\!\left(z,q\right)+\sum_{n=1% }^{\infty}\left(\alpha_{n}\mathop{\mathrm{ce}_{n}\/}\nolimits\!\left(z,q\right% )+\beta_{n}\mathop{\mathrm{se}_{n}\/}\nolimits\!\left(z,q\right)\right),$ Defines: $f(z)$: function (locally) Symbols: $\mathop{\mathrm{ce}_{\NVar{n}}\/}\nolimits\!\left(\NVar{z},\NVar{q}\right)$: Mathieu function, $\mathop{\mathrm{se}_{\NVar{n}}\/}\nolimits\!\left(\NVar{z},\NVar{q}\right)$: Mathieu function, $q=h^{2}$: parameter, $n$: integer, $z$: complex variable, $\alpha_{n}$ and $\beta_{n}$ Referenced by: §28.11 Permalink: http://dlmf.nist.gov/28.11.E1 Encodings: TeX, pMML, png See also: Annotations for 28.11

where

 28.11.2 $\displaystyle\alpha_{n}$ $\displaystyle=\frac{1}{\pi}\int_{0}^{2\pi}f(x)\mathop{\mathrm{ce}_{n}\/}% \nolimits\!\left(x,q\right)\mathrm{d}x,$ $\displaystyle\beta_{n}$ $\displaystyle=\frac{1}{\pi}\int_{0}^{2\pi}f(x)\mathop{\mathrm{se}_{n}\/}% \nolimits\!\left(x,q\right)\mathrm{d}x.$

The series (28.11.1) converges absolutely and uniformly on any compact subset of the strip $S$. See Meixner and Schäfke (1954, §2.28), and for expansions in the case of the exceptional values of $q$ see Meixner et al. (1980, p. 33).

## Examples

With the notation of §28.4,

 28.11.3 $1=2\sum_{n=0}^{\infty}A_{0}^{2n}(q)\mathop{\mathrm{ce}_{2n}\/}\nolimits\!\left% (z,q\right),$
 28.11.4 $\mathop{\cos\/}\nolimits 2mz=\sum_{n=0}^{\infty}A_{2m}^{2n}(q)\mathop{\mathrm{% ce}_{2n}\/}\nolimits\!\left(z,q\right),$ $m\neq 0$,
 28.11.5 $\displaystyle\mathop{\cos\/}\nolimits(2m+1)z$ $\displaystyle=\sum_{n=0}^{\infty}A_{2m+1}^{2n+1}(q)\mathop{\mathrm{ce}_{2n+1}% \/}\nolimits\!\left(z,q\right),$ 28.11.6 $\displaystyle\mathop{\sin\/}\nolimits(2m+1)z$ $\displaystyle=\sum_{n=0}^{\infty}B_{2m+1}^{2n+1}(q)\mathop{\mathrm{se}_{2n+1}% \/}\nolimits\!\left(z,q\right),$ 28.11.7 $\displaystyle\mathop{\sin\/}\nolimits(2m+2)z$ $\displaystyle=\sum_{n=0}^{\infty}B_{2m+2}^{2n+2}(q)\mathop{\mathrm{se}_{2n+2}% \/}\nolimits\!\left(z,q\right).$


text:
# 3.11 Molecular orbital theory  (Page 6/26)

 Page 6 / 26
$\text{bond order}=\phantom{\rule{0.2em}{0ex}}\frac{\left(\text{number of bonding electrons}\right)\phantom{\rule{0.2em}{0ex}}-\phantom{\rule{0.2em}{0ex}}\left(\text{number of antibonding electrons}\right)}{2}$

The order of a covalent bond is a guide to its strength; a bond between two given atoms becomes stronger as the bond order increases ( [link] ). If the distribution of electrons in the molecular orbitals between two atoms is such that the resulting bond would have a bond order of zero, a stable bond does not form. We next look at some specific examples of MO diagrams and bond orders.

## Bonding in diatomic molecules

A dihydrogen molecule (H 2 ) forms from two hydrogen atoms. When the atomic orbitals of the two atoms combine, the electrons occupy the molecular orbital of lowest energy, the σ 1 s bonding orbital. A dihydrogen molecule, H 2 , readily forms because the energy of a H 2 molecule is lower than that of two H atoms. The σ 1 s orbital that contains both electrons is lower in energy than either of the two 1 s atomic orbitals.

A molecular orbital can hold two electrons, so both electrons in the H 2 molecule are in the σ 1 s bonding orbital; the electron configuration is ${\left({\text{σ}}_{1s}\right)}^{2}.$ We represent this configuration by a molecular orbital energy diagram ( [link] ) in which a single upward arrow indicates one electron in an orbital, and two (upward and downward) arrows indicate two electrons of opposite spin.

A dihydrogen molecule contains two bonding electrons and no antibonding electrons so we have

${\text{bond order in H}}_{2}=\phantom{\rule{0.2em}{0ex}}\frac{\left(2-0\right)}{2}\phantom{\rule{0.2em}{0ex}}=1$

Because the bond order for the H–H bond is equal to 1, the bond is a single bond.

A helium atom has two electrons, both of which are in its 1 s orbital. Two helium atoms do not combine to form a dihelium molecule, He 2 , with four electrons, because the stabilizing effect of the two electrons in the lower-energy bonding orbital would be offset by the destabilizing effect of the two electrons in the higher-energy antibonding molecular orbital. We would write the hypothetical electron configuration of He 2 as ${\left({\text{σ}}_{1s}\right)}^{2}{\left({\text{σ}}_{1s}^{*}\right)}^{2}$ as in [link] . The net energy change would be zero, so there is no driving force for helium atoms to form the diatomic molecule. In fact, helium exists as discrete atoms rather than as diatomic molecules. The bond order in a hypothetical dihelium molecule would be zero.

${\text{bond order in He}}_{2}=\phantom{\rule{0.2em}{0ex}}\frac{\left(2-2\right)}{2}\phantom{\rule{0.2em}{0ex}}=0$

A bond order of zero indicates that no bond is formed between two atoms.

## The diatomic molecules of the second period

Eight possible homonuclear diatomic molecules might be formed by the atoms of the second period of the periodic table: Li 2 , Be 2 , B 2 , C 2 , N 2 , O 2 , F 2 , and Ne 2 . However, we can predict that the Be 2 molecule and the Ne 2 molecule would not be stable. We can see this by a consideration of the molecular electron configurations ( [link] ).

We predict valence molecular orbital electron configurations just as we predict electron configurations of atoms. Valence electrons are assigned to valence molecular orbitals with the lowest possible energies. Consistent with Hund’s rule, whenever there are two or more degenerate molecular orbitals, electrons fill each orbital of that type singly before any pairing of electrons takes place.

#### Questions & Answers

I only see partial conversation and what's the question here!
what about nanotechnology for water purification
please someone correct me if I'm wrong but I think one can use nanoparticles, specially silver nanoparticles for water treatment.
Damian
what is the stm
is there industrial application of fullrenes. What is the method to prepare fullrene on large scale.?
Rafiq
industrial application...? mmm I think on the medical side as drug carrier, but you should go deeper on your research, I may be wrong
Damian
How we are making nano material?
what is a peer
What is meant by 'nano scale'?
What is STMs full form?
LITNING
scanning tunneling microscope
Sahil
how nano science is used for hydrophobicity
Santosh
Do u think that Graphene and Fullrene fiber can be used to make Air Plane body structure the lightest and strongest. Rafiq
Rafiq
what is differents between GO and RGO?
Mahi
what is simplest way to understand the applications of nano robots used to detect the cancer affected cell of human body.? How this robot is carried to required site of body cell.? what will be the carrier material and how can be detected that correct delivery of drug is done Rafiq
Rafiq
what is Nano technology ?
write examples of Nano molecule?
Bob
The nanotechnology is as new science, to scale nanometric
brayan
nanotechnology is the study, desing, synthesis, manipulation and application of materials and functional systems through control of matter at nanoscale
Damian
Is there any normative that regulates the use of silver nanoparticles?
what king of growth are you checking .?
Renato
What fields keep nano created devices from performing or assimulating ? Magnetic fields ? Are do they assimilate ?
why we need to study biomolecules, molecular biology in nanotechnology?
?
Kyle
yes I'm doing my masters in nanotechnology, we are being studying all these domains as well..
why?
what school?
Kyle
biomolecules are e building blocks of every organics and inorganic materials.
Joe
anyone know any internet site where one can find nanotechnology papers?
research.net
kanaga
sciencedirect big data base
Ernesto
Introduction about quantum dots in nanotechnology
what does nano mean?
nano basically means 10^(-9). nanometer is a unit to measure length.
Bharti
do you think it's worthwhile in the long term to study the effects and possibilities of nanotechnology on viral treatment?
absolutely yes
Daniel
how to know photocatalytic properties of tio2 nanoparticles...what to do now
it is a goid question and i want to know the answer as well
Maciej
characteristics of micro business
Abigail
for teaching engĺish at school how nano technology help us
Anassong
How can I make nanorobot?
Lily
how did you get the value of 2000N.What calculations are needed to arrive at it
Privacy Information Security Software Version 1.1a
Good
how do you find theWhat are the wavelengths and energies per photon of two lines
The eyes of some reptiles are sensitive to 850 nm light. If the minimum energy to trigger the receptor at this wavelength is 3.15 x 10-14 J, what is the minimum number of 850 nm photons that must hit the receptor in order for it to be triggered?
A teaspoon of the carbohydrate sucrose contains 16 calories, what is the mass of one teaspoo of sucrose if the average number of calories for carbohydrate is 4.1 calories/g?
4. On the basis of dipole moments and/or hydrogen bonding, explain in a qualitative way the differences in the boiling points of acetone (56.2 °C) and 1-propanol (97.4 °C), which have similar molar masses
Calculate the bond order for an ion with this configuration: (?2s)2(??2s)2(?2px)2(?2py,?2pz)4(??2py,??2pz)3
Which of the following will increase the percent of HF that is converted to the fluoride ion in water? (a) addition of NaOH (b) addition of HCl (c) addition of NaF