Repository: ecmwf/thermofeel (language: python)

Issue:
Reassignment throws error in calculate_cos_solar_zenith_angle_integrated
## The Bug

While attempting to use calculate_cos_solar_zenith_angle_integrated via xr.apply_ufunc as in: 
```python
integralCzaEaHour = [xr.apply_ufunc(thermofeel.calculate_cos_solar_zenith_angle_integrated, ds.lat, ds.lon, time.dt.year, 
    time.dt.month, time.dt.day, time.dt.hour, 0, 1) for time in ds.time]
```
I get the following error:
`ValueError: non-broadcastable output operand with shape (571,1) doesn't match the broadcast shape (571,1440)`

## The Solution

An [answer](https://stackoverflow.com/a/47493938/8891677) to a related stack overflow question suggests eliminating the `+=` operator, as in changing `a += b` to `a = a + b`. Sure enough, when I change [line 316 in thermofeel.py](https://github.com/ecmwf-projects/thermofeel/blob/720911162ba1c614a20d6fd3c87dea9f554c5fd0/thermofeel/thermofeel.py#L316) from `integral += w[n] * calculate_cos_solar_zenith_angle(` to `integral = integral + w[n] * calculate_cos_solar_zenith_angle(`, everything works as expected. The stack overflow answer includes an explanation as to why.

Better explanation of T0
This is related to issue #34 which was implemented in 2.0.0 but with slightly different values for T0, the triple point of water. 
fix: fix mistake in `normal_effective_temperature()`
I think you missed the brackets

https://link.springer.com/article/10.1007/s00484-011-0453-2


Discussion / hints:
@jarataraj many thanks for pointing out this issue. 
I think it makes sense to apply the fix. Would you be able to provide a pull request with solution?
@tlmquintino Will do! Is it ok for me to change all instances of `a += b` and similar to `a = a + b` and similar? Or should I just apply the fix to `calculate_cos_solar_zenith_angle_integrated`? I know I got the error for two different thermofeel functions but I can't remember what the other was and haven't yet tested all of them. 
@jarataraj sorry for the late reply, I've been traveling. 
Yes please do go ahead and I think it makes sense doing that substitution everywhere.



Fix patch:
diff --git a/ChangeLog.rst b/ChangeLog.rst
index ba334ca..5fff7e2 100644
--- a/ChangeLog.rst
+++ b/ChangeLog.rst
@@ -1,6 +1,12 @@
 ChangeLog
 =========
 
+Version 2.1.0
+-------------
+
+ * Fix computaiton of normal effective temperature
+ * Fix documentation links to new repository URL
+ 
 Version 2.0.0
 -------------
 
diff --git a/README.rst b/README.rst
index 42dcbcd..b3d7d61 100644
--- a/README.rst
+++ b/README.rst
@@ -1,8 +1,8 @@
-.. image:: https://raw.githubusercontent.com/ecmwf-projects/thermofeel/master/thermofeel.png
+.. image:: https://raw.githubusercontent.com/ecmwf/thermofeel/master/thermofeel.png
   :width: 600
   :alt: thermofeel logo
 
-|license| |tag_release| |docs|
+|license| |tag_release| |docs| |ci|
 
 **thermofeel** (pronounced *thermo-feel*)
 
@@ -36,11 +36,6 @@ Install with::
 
     $ pip install thermofeel
 
-Testing
-=======
-
-|ci| |codecov|
-
 System dependencies
 ===================
 
@@ -68,24 +63,24 @@ The main changes are:
 
 Please consult ChangeLog_ for more details.
 
-.. _ChangeLog: https://github.com/ecmwf-projects/thermofeel/blob/master/ChangeLog.rst
+.. _ChangeLog: https://github.com/ecmwf/thermofeel/blob/master/ChangeLog.rst
 
 
 Contributing
 ============
 
-The main repository is hosted on `GitHub <https://github.com/ecmwf-projects/thermofeel>`_. Testing, bug reports and contributions are highly welcomed and appreciated.
+The main repository is hosted on `GitHub <https://github.com/ecmwf/thermofeel>`_. Testing, bug reports and contributions are highly welcomed and appreciated.
 
 Please see the Contributing_ document for the best way to help.
 
-.. _Contributing: https://github.com/ecmwf-projects/thermofeel/blob/master/CONTRIBUTING.rst
+.. _Contributing: https://github.com/ecmwf/thermofeel/blob/master/CONTRIBUTING.rst
 
 Current developers:
 
 - Claudia Di Napoli - `ECMWF <https://ecmwf.int>`_
 - Tiago Quintino - `ECMWF <https://ecmwf.int>`_
 
-See also the `contributors <https://github.com/ecmwf-projects/thermofeel/contributors>`_ for a more complete list.
+See also the `contributors <https://github.com/ecmwf/thermofeel/contributors>`_ for a more complete list.
 
 License
 =======
@@ -130,16 +125,16 @@ Acknowledgements
 Past and current funding and support for **thermofeel** is listed in the adjoning Acknowledgements_
 
 
-.. _Acknowledgements: https://github.com/ecmwf-projects/thermofeel/blob/master/ACKNOWLEDGEMENTS.rst
+.. _Acknowledgements: https://github.com/ecmwf/thermofeel/blob/master/ACKNOWLEDGEMENTS.rst
 
 
-.. |last_commit| image:: https://img.shields.io/github/last-commit/ecmwf-projects/thermofeel
-    :target: https://github.com/ecmwf-projects/thermofeel
+.. |last_commit| image:: https://img.shields.io/github/last-commit/ecmwf/thermofeel
+    :target: https://github.com/ecmwf/thermofeel
 
-.. |commits_since_release| image:: https://img.shields.io/github/commits-since/ecmwf-projects/thermofeel/latest?sort=semver
-    :target: https://github.com/ecmwf-projects/thermofeel
+.. |commits_since_release| image:: https://img.shields.io/github/commits-since/ecmwf/thermofeel/latest?sort=semver
+    :target: https://github.com/ecmwf/thermofeel
 
-.. |license| image:: https://img.shields.io/github/license/ecmwf-projects/thermofeel
+.. |license| image:: https://img.shields.io/github/license/ecmwf/thermofeel
     :target: https://www.apache.org/licenses/LICENSE-2.0.html
 
 .. |pypi_release| image:: https://img.shields.io/pypi/v/thermofeel?color=green
@@ -148,20 +143,20 @@ Past and current funding and support for **thermofeel** is listed in the adjonin
 .. |pypi_status| image:: https://img.shields.io/pypi/status/thermofeel
     :target: https://pypi.org/project/thermofeel
 
-.. |tag_release| image:: https://img.shields.io/github/v/release/ecmwf-projects/thermofeel?sort=semver
-    :target: https://github.com/ecmwf-projects/thermofeel
+.. |tag_release| image:: https://img.shields.io/github/v/release/ecmwf/thermofeel?sort=semver
+    :target: https://github.com/ecmwf/thermofeel
 
-.. |codecov| image:: https://codecov.io/gh/ecmwf-projects/thermofeel/branch/master/graph/badge.svg
-  :target: https://codecov.io/gh/ecmwf-projects/thermofeel
+.. |codecov| image:: https://codecov.io/gh/ecmwf/thermofeel/branch/master/graph/badge.svg
+  :target: https://codecov.io/gh/ecmwf/thermofeel
 
-.. |ci| image:: https://img.shields.io/github/workflow/status/ecmwf-projects/thermofeel/ci
-  :target: https://github.com/ecmwf-projects/thermofeel/actions
+.. |ci| image:: https://img.shields.io/github/actions/workflow/status/ecmwf/thermofeel/ci.yml
+  :target: https://github.com/ecmwf/thermofeel/actions
 
 .. |pypi_downloads| image:: https://img.shields.io/pypi/dm/thermofeel
   :target: https://pypi.org/project/thermofeel
 
-.. |code_size| image:: https://img.shields.io/github/languages/code-size/ecmwf-projects/thermofeel?color=green
-  :target: https://github.com/ecmwf-projects/thermofeel
+.. |code_size| image:: https://img.shields.io/github/languages/code-size/ecmwf/thermofeel?color=green
+  :target: https://github.com/ecmwf/thermofeel
   
 .. |docs| image:: https://readthedocs.org/projects/thermofeel/badge/?version=latest
   :target: https://thermofeel.readthedocs.io/en/latest/?badge=latest
diff --git a/thermofeel/__init__.py b/thermofeel/__init__.py
index cb620eb..6ecce3a 100644
--- a/thermofeel/__init__.py
+++ b/thermofeel/__init__.py
@@ -8,4 +8,4 @@
 
 from .thermofeel import *  # noqa
 
-__version__ = "2.0.0"
+__version__ = "2.1.0"
diff --git a/thermofeel/thermofeel.py b/thermofeel/thermofeel.py
index 90854c2..7a976b7 100644
--- a/thermofeel/thermofeel.py
+++ b/thermofeel/thermofeel.py
@@ -85,7 +85,7 @@ def calculate_saturation_vapour_pressure(t2_k):
     ]
     ess = g[7] * np.log(t2_k)
     for i in range(7):
-        ess += g[i] * np.power(t2_k, (i - 2))
+        ess = ess + g[i] * np.power(t2_k, (i - 2))
 
     ess = np.exp(ess) * 0.01  # hPa
 
@@ -98,14 +98,15 @@ def calculate_saturation_vapour_pressure_multiphase(t2_k, phase):
         :param t2_k: (float array) 2m temperature [K]
         :param phase: 0 over liquid water and 1 over ice
         returns pressure of water vapor over a surface of liquid water or ice [hPa] == [mBar]
-    Reference: ECMWF IFS Documentation CY45R1 - Part IV : Physical processes (2018)
+    Reference: ECMWF IFS Documentation CY45R1 - Part IV : Physical processes (2018) pp. 116
     https://doi.org/10.21957/4whwo8jw0
+    https://metview.readthedocs.io/en/latest/api/functions/saturation_vapour_pressure.html
     """
-
+    T0 = 273.16  # triple point of water 273.16 K (0.01 °C) at 611.73 Pa
     es = np.zeros_like(t2_k)
-    y = (t2_k - 273.16) / (t2_k - 32.19)  # over liquid water
+    y = (t2_k - T0) / (t2_k - 32.19)  # over liquid water
     es[phase == 0] = 6.1121 * np.exp(17.502 * y[phase == 0])
-    y = (t2_k - 273.16) / (t2_k + 0.7)  # over ice
+    y = (t2_k - T0) / (t2_k + 0.7)  # over ice
     es[phase == 1] = 6.1121 * np.exp(22.587 * y[phase == 1])
 
     return es
@@ -677,7 +678,7 @@ def calculate_normal_effective_temperature(t2_k, va, rh):
     ditermeq = 1 / (1.76 + 1.4 * v**0.75)
     net = (
         37
-        - (37 - t2_k / (0.68 - 0.0014 * rh + ditermeq))
+        - ((37 - t2_k) / (0.68 - 0.0014 * rh + ditermeq))
         - 0.29 * t2_k * (1 - 0.01 * rh)
     )
     net_k = celsius_to_kelvin(net)
@@ -710,11 +711,11 @@ def calculate_wind_chill(t2_k, va):
         :param t2_k: (float array) 2m Temperature [K]
         :param va: (float array) wind speed at 10 meters [m/s]
         returns wind chill [K]
-        Temperature must be between -50°C and 5°C; wind speed must be between 5km/h and 80km/h
-        Wind chill from input values outside those ranges are not to be considered valid
+        Computation is only valid for temperatures between -50°C and 5°C and wind speeds between 5km/h and 80km/h.
+        For input values outside those ranges, computed results not be considered valid.
     Reference: Blazejczyk et al. (2012)
     https://doi.org/10.1007/s00484-011-0453-2
-    See also: http://www.ec.gc.ca/meteo-weather/default.asp?lang=n&n=5FBF816A-1#wc6
+    See also: https://web.archive.org/web/20130627223738/http://climate.weatheroffice.gc.ca/prods_servs/normals_documentation_e.html  # noqa
     """
     t2_c = kelvin_to_celsius(t2_k)  # kelvin_to_celsius(tk)
     v = va * 3.6  # convert to kilometers per hour


Test patch:
diff --git a/tests/net.csv b/tests/net.csv
index a42e4e7..7ab9b8b 100644
--- a/tests/net.csv
+++ b/tests/net.csv
@@ -1,50 +1,50 @@
-3.033889702205782442e+02
-2.928502091976476436e+02
-2.770285138080615752e+02
-2.770490005238265780e+02
-2.770928248276744625e+02
-2.771638191543268022e+02
-2.772180430167822465e+02
-2.773481405053552180e+02
-2.774187263366436014e+02
-2.774623803824222819e+02
-2.775973885081999697e+02
-2.776094978681781527e+02
-2.775812180858250713e+02
-2.775939616382239024e+02
-2.775908341477433510e+02
-2.775773393811281267e+02
-2.775198416488409521e+02
-2.774920248106813006e+02
-2.774635394202203997e+02
-2.773664551403895189e+02
-2.773095070329615623e+02
-2.772699512899097840e+02
-2.772017131727243964e+02
-2.771203280402834821e+02
-2.769705081825984507e+02
-2.769574995185245712e+02
-2.937299618046274645e+02
-2.935924042492049466e+02
-2.933120777011570226e+02
-2.932935807975503053e+02
-2.933526871192943872e+02
-2.936544161860668964e+02
-2.938921675892892722e+02
-2.941795452139757003e+02
-2.947574518545506521e+02
-2.950281071095799916e+02
-2.950420555947761159e+02
-2.948682023918578352e+02
-2.946889368890751939e+02
-2.945812316730967950e+02
-2.945970172400211027e+02
-2.945403954938864217e+02
-2.947040000154167956e+02
-2.947692211015685757e+02
-2.949256961037862652e+02
-2.951105589804234341e+02
-2.949321900709687725e+02
-2.947569560678213065e+02
-2.942588330787240238e+02
-2.942205457588216859e+02
+3.010065622860677763e+02
+2.978004713493879763e+02
+2.762241215897753364e+02
+2.762341286739048769e+02
+2.763012956155728830e+02
+2.764059939851455852e+02
+2.764790654938740317e+02
+2.766273813534896249e+02
+2.766893014126333696e+02
+2.766966993523136011e+02
+2.767887884697527738e+02
+2.767496737064529952e+02
+2.769510933062428535e+02
+2.768837969370093788e+02
+2.768171117155841898e+02
+2.767348212481284122e+02
+2.766230919517424240e+02
+2.765545420051105907e+02
+2.765267374915458731e+02
+2.764466352440045398e+02
+2.764088648567151267e+02
+2.763790656581484200e+02
+2.763100466759992742e+02
+2.762046876422498372e+02
+2.762623057754365163e+02
+2.761689174526972579e+02
+2.911413948240829654e+02
+2.910768702461481325e+02
+2.907986242602605671e+02
+2.907857526147024601e+02
+2.907969231665560415e+02
+2.911296388794274321e+02
+2.914604298204376391e+02
+2.918874048860148491e+02
+2.928703177480205682e+02
+2.936264195426511492e+02
+2.948825390578576844e+02
+2.951751610936696579e+02
+2.953622575333719738e+02
+2.954328487477691283e+02
+2.953756039592274760e+02
+2.949602578417895415e+02
+2.946793060090293466e+02
+2.939725759663496092e+02
+2.934531765165319825e+02
+2.931004552660506874e+02
+2.923946543653016761e+02
+2.918684112352694342e+02
+2.913048255487498182e+02
+2.910969750368326459e+02
diff --git a/tests/test_scalars.py b/tests/test_scalars.py
index e13f716..a1f5eaa 100644
--- a/tests/test_scalars.py
+++ b/tests/test_scalars.py
@@ -144,7 +144,7 @@ def test_normal_effective_temperature(self):
         rh = np.array([80])
         net = np.array([tmf.calculate_normal_effective_temperature(t2_k, va, rh)])
         # print(f"net {net}")
-        assert net == pytest.approx(314.7102642987715, abs=1e-6)
+        assert net == pytest.approx(304.13650125, abs=1e-6)
 
     def test_apparent_temperature(self):
         t2_k = np.array([tmf.celsius_to_kelvin(25.0)])
@@ -154,11 +154,26 @@ def test_apparent_temperature(self):
         # print(f"at {at}")
         assert at == pytest.approx(299.86678322384626, abs=1e-6)
 
-    def test_windchill(self):
+    def test_wind_chill(self):
         t2_k = np.array([270])
         va = np.array([10])
-        wc = np.array([tmf.calculate_wind_chill(t2_k, va)])
-        assert wc == pytest.approx(261.92338925380074, abs=1e-6)
+        wc_k = np.array([tmf.calculate_wind_chill(t2_k, va)])
+        assert wc_k == pytest.approx(261.92338925380074, abs=1e-6)
+        # reference result from from wikipedia article https://en.wikipedia.org/wiki/Wind_chill
+        t2_k = np.array([tmf.celsius_to_kelvin(-20)])  # -20C to K
+        va = np.array([5 / 3.6])  # 5 km/h to m/s
+        wc_k = np.array([tmf.calculate_wind_chill(t2_k, va)])
+        wc_c = tmf.kelvin_to_celsius(wc_k)
+        assert wc_c == pytest.approx(
+            -24.27850328, abs=1e-6
+        )  # around ~ -24C for wind chill (not exact)
+        va = np.array([30 / 3.6])  # 30 km/h to m/s
+        wc_k = np.array([tmf.calculate_wind_chill(t2_k, va)])
+        wc_c = tmf.kelvin_to_celsius(wc_k)
+        print(f"wc_c {wc_c}")
+        assert wc_c == pytest.approx(
+            -32.56804448, abs=1e-6
+        )  # around ~ -33C for wind chill (not exact)
 
     def test_heat_index_simplified(self):
         t2_k = np.array([tmf.celsius_to_kelvin(21.0)])