Repository: fatiando/boule

Issue:
Check if normal gravity is computed on internal point
**Description of the desired feature**

The normal gravity formulas used on `boule.Ellipsoid.normal_gravity()` and `boule.Sphere.normal_gravity()` are valid only for points outside the ellipsoid.
Would be nice to add a warning in case the given point falls inside it.
To do so we just need to check if the `height` is positive or equal to zero.

**Are you willing to help implement and maintain this feature?** Yes, but would be glad to help anyone who wants to implement it.

<!--
Every feature we add is code that we will have to maintain and keep updated.
This takes a lot of effort. If you are willing to be involved in the project
and help maintain your feature, it will make it easier for us to accept it.
-->


Discussion:
MGomezN: Hi! I am back. I just wanted to make a pull request for my last commit, but I couldn't erase the past commits. This time I tried to follow more carefully the Contributing Guidelines. I ran checks format,  check,  lint and test. All of them appear successful. However, I am dubious about the solution I implemented using the norm. Maybe is a condition that is always true by definition of "norm"... so, I'll wait for your review. Thank you!
santisoler: HI @MGomezN! Thanks again for your contributions.

> I just wanted to make a pull request for my last commit, but I couldn't erase the past commits.

I'm not sure if I'm following you here. Have you tried to commit to your `internal_point` branch? In git we don't erase previous commits (unless it's strictly necessary) because that could potentially create a lot of troubles. Instead you can always add a new commit making changes that fixes the problems introduced by any previous commit. This not only makes version control easier, but also helps us to keep track of the changes and decisions we made in the past.

> However, I am dubious about the solution I implemented using the norm. Maybe is a condition that is always true by definition of "norm"... so, I'll wait for your review.

I'm not entirely sure why are you checking if the norm of the height is greater than zero. Raising the warning if `height < 0` should be enough, right?

On #66 you also added the same check for the `Sphere` class. Would you like to add it here if you want to continue the development of this feature on this PR? 
MGomezN: > HI @MGomezN! Thanks again for your contributions.
> 
 Thank you @santisoler 
>
> I'm not sure if I'm following you here. Have you tried to commit to your `internal_point` branch? In git, we don't erase previous commits (unless it's strictly necessary) because that could potentially create a lot of troubles. Instead you can always add a new commit making changes that fixes the problems introduced by any previous commit. This not only makes version control easier, but also helps us to keep track of the changes and decisions we made in the past.
> 
No. I didn't try to commit to my internal_point branch. I worked with internal_point branch almost six months ago, and by that time, I had a lot of trouble with the environment, and also I made a reset of my PC. Therefore It seemed easier to start everything from zero last week when I resume my work on this issue. So, I closed my pull request related to the internal_point branch, and actually, I erased my fork of the Boule repository, and I started everything from the beginning: clone, fork, Install environment, and so on. I incorporated comments from Leonardo in internal_point branch, and as a result, I created commit 9b71364.
>
But… I forgot to run make format, make check, make linnet, so I did it. I discovered one error, I repaired that error, and then I created commit bd0dcd7.
>
> I'm not entirely sure why are you checking if the norm of the height is greater than zero. Raising the warning if `height < 0` should be enough, right?
>
Again, I forgot to ran the make test, I did it, and there was a new error saying height <0 is not ok because height could be an array, not just a scalar, so I tried to use the norm of the array “height,” and that was mi final commit  a0a2239
>
So, I wanted to erase commits 9b71364 and bd0dcd7 because I know (now) they were incomplete. The most complete is commit a0a2239.
Jeje It seems I still don't get a good understanding of Git, but for me, it is unnecessary to leave 9b71364 and  bd0dcd7 in the track of changes. That is why I said I wanted to delete them.
> 
> On #66 you also added the same check for the `Sphere` class. Would you like to add it here if you want to continue the development of this feature on this PR?

Ups, see? Not too much focus in my head. Yes, I want to change the Sphere class also.
santisoler: > No. I didn't try to commit to my internal_point branch. I worked with internal_point branch almost six months ago, and by that time, I had a lot of trouble with the environment, and also I made a reset of my PC. Therefore It seemed easier to start everything from zero last week when I resume my work on this issue. So, I closed my pull request related to the internal_point branch, and actually, I erased my fork of the Boule repository, and I started everything from the beginning: clone, fork, Install environment, and so on. I incorporated comments from Leonardo made on the internal_point branch, and as a result, I created commit [9b71364](https://github.com/fatiando/boule/commit/9b713645adcc88b397b3919718f6f58a16fa067f).
> But… I forgot to run make format, make check, make linnet, so I did it. I discovered one error, I repaired that error, and then I created commit [bd0dcd7](https://github.com/fatiando/boule/commit/bd0dcd7358a23a96dfca6006bff0eec7d828d46a).

No big deal then! Just wanted to know if you were experiencing technical issues or if you willingly started again from scratch.

> Again, I forgot to ran the make test, I did it, and there was a new error saying height <0 is not ok because height could be an array, not just a scalar, so I tried to use the norm of the array “height,” and that was mi final commit [a0a2239](https://github.com/fatiando/boule/commit/a0a22394281e99f2aecb44e158c24f527bf55dc6)

You are absolutely right! We cannot use `height < 0` because `height` can be an array and it would return a boolean array, so the statement `height < 0` cannot be evaluated as a single bool. Nevertheless, using `np.linalg.norm` would create always positive values, so it won't work for our purposes. What we need to check if *any* value of `height` is negative, so give `np.any()` a try.

> So, I wanted to erase commits [9b71364](https://github.com/fatiando/boule/commit/9b713645adcc88b397b3919718f6f58a16fa067f) and [bd0dcd7](https://github.com/fatiando/boule/commit/bd0dcd7358a23a96dfca6006bff0eec7d828d46a) because I know (now) they were incomplete. The most complete is commit [a0a2239](https://github.com/fatiando/boule/commit/a0a22394281e99f2aecb44e158c24f527bf55dc6).
> Jeje It seems I still don't get a good understanding of Git, but for me, it is unnecessary to leave [9b71364](https://github.com/fatiando/boule/commit/9b713645adcc88b397b3919718f6f58a16fa067f) and [bd0dcd7](https://github.com/fatiando/boule/commit/bd0dcd7358a23a96dfca6006bff0eec7d828d46a) in the track of changes. That is why I said I wanted to delete them.

We virtually never *delete* a commit, as I said before, the best practice is to add a new commit that fixes the mistakes of the previous ones. Think of git as a journey: you know where you started and you know where you are right now, but you cannot erase all the points in the middle, because they took you where you are right now.
That being said, when we merge this branch into `master` we do a Squash and Merge: we squash every commit from this PR into a single one. So feel free to add as many commits as you want, they will not end up in the history of Boule's `master` branch.
Using Squash and Merge generates a way cleaner log history on `master` and makes the development of new features much easier by allowing us to break stuff on any other branch, we just need to be very careful what we merge to `master` and that's why we have all the CIs automatically checking stuff for us.

> Ups, see? Not too much focus in my head. Yes, I want to change the Sphere class also.

No problem! Feel free to finish the warning on `Ellipsoid` and once you got it right, you can start adding it to `Sphere`.

MGomezN: > You are absolutely right! We cannot use `height < 0` because `height` can be an array and it would return a boolean array, so the statement `height < 0` cannot be evaluated as a single bool. Nevertheless, using `np.linalg.norm` would create always positive values, so it won't work for our purposes. What we need to check if _any_ value of `height` is negative, so give `np.any()` a try.
>
Hello @santisoler. I tried your suggestion in commit 8b48925
> We virtually never _delete_ a commit, as I said before, the best practice is to add a new commit that fixes the mistakes of the previous ones. Think of git as a journey: you know where you started and you know where you are right now, but you cannot erase all the points in the middle, because they took you where you are right now.
> That being said, when we merge this branch into `master` we do a Squash and Merge: we squash every commit from this PR into a single one. So feel free to add as many commits as you want, they will not end up in the history of Boule's `master` branch.
> Using Squash and Merge generates a way cleaner log history on `master` and makes the development of new features much easier by allowing us to break stuff on any other branch, we just need to be very careful what we merge to `master` and that's why we have all the CIs automatically checking stuff for us.
>
I understand.  Looking forward to squash and merge all my mistakes jeje
>
> No problem! Feel free to finish the warning on `Ellipsoid` and once you got it right, you can start adding it to `Sphere.`
>
I think it works now. I only feel insecure about the place where I wrote the warning (under de def of normal_gravity). I don't know if it must be written before, for example, after line 97, where are those bunch of `checks` e.g. `@flattening.validator`. 

santisoler: > I understand. Looking forward to squash and merge all my mistakes jeje

😂

> I think it works now. I only feel insecure about the place where I wrote the warning (under de def of normal_gravity). I don't know if it must be written before, for example, after line 97, where are those bunch of `checks` e.g. `@flattening.validator`.

The warning is coded in the perfect spot. The requirement of `height` being positive or zero only applies to the `normal_gravity` method, so there's no point in creating a separate private method (like `_check_flattening`) for raising this particular warning: no other method will use it, and because the warning only takes ~3 lines, we are not cluttering the code.

Would you like to start writing a test function that checks if the warning is raised after passing a negative height? 
If so, please create a new function in `tests/test_ellipsoid.py`. You can catch warnings in test functions with `warnings.catch_warnings(record=True)`, like this:
https://github.com/fatiando/boule/blob/145384d87a68569dc4ede3acd2a39bbffe16c6b5/boule/tests/test_ellipsoid.py#L59-L67

If you need any help, please feel free to ask!
MGomezN: > Would you like to start writing a test function that checks if the warning is raised after passing a negative height?
> If so, please create a new function in `tests/test_ellipsoid.py`. You can catch warnings in test functions with `warnings.catch_warnings(record=True)`, like this:
> https://github.com/fatiando/boule/blob/145384d87a68569dc4ede3acd2a39bbffe16c6b5/boule/tests/test_ellipsoid.py#L59-L67
> 
Thank you @santisoler . I added the test and the warning for boule.Sphere.normal_gravity() in commit 1f6f555. Fingers cross!  

leouieda: And a big thank you to @santisoler for providing brilliant guidance here (as always) 🎉 

Patch:
diff --git a/boule/ellipsoid.py b/boule/ellipsoid.py
index b75498b..a71cf6a 100644
--- a/boule/ellipsoid.py
+++ b/boule/ellipsoid.py
@@ -453,6 +453,13 @@ class Ellipsoid:
             The normal gravity in mGal.
 
         """
+        # Warn if height is negative
+        if np.any(height < 0):
+            warn(
+                "Formulas used are valid for points outside the ellipsoid."
+                "Height must be greater than or equal to zero."
+            )
+
         sinlat = np.sin(np.deg2rad(latitude))
         coslat = np.sqrt(1 - sinlat ** 2)
         # The terms below follow the variable names from Li and Goetze (2001)
diff --git a/boule/sphere.py b/boule/sphere.py
index 15e8ea7..895cadd 100644
--- a/boule/sphere.py
+++ b/boule/sphere.py
@@ -193,6 +193,12 @@ class Sphere(Ellipsoid):
             The normal gravity in mGal.
 
         """
+        # Warn if height is negative
+        if np.any(height < 0):
+            warn(
+                "Formulas used are valid for points outside the ellipsoid."
+                "Height must be greater than or equal to zero."
+            )
         radial_distance = self.radius + height
         gravity_acceleration = self.geocentric_grav_const / (radial_distance) ** 2
         gamma = np.sqrt(


Test patch:
diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py
index b509f43..e96dabf 100644
--- a/boule/tests/test_ellipsoid.py
+++ b/boule/tests/test_ellipsoid.py
@@ -325,3 +325,14 @@ def test_normal_gravity_against_somigliana(ellipsoid):
             ellipsoid.normal_gravity(latitude, height=0),
             normal_gravity_surface(latitude, ellipsoid),
         )
+
+
+@pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)
+def test_normal_gravity_computed_on_internal_point(ellipsoid):
+    """
+    Check if warn is raised if height is negative
+    """
+    latitude = np.linspace(-90, 90, 100)
+    with warnings.catch_warnings(record=True) as warn:
+        ellipsoid.normal_gravity(latitude, height=-10)
+        assert len(warn) >= 1