TEOS-10 / GibbsSeaWater.jl

Gibbs-SeaWater (GSW) Oceanographic Toolbox in Julia
http://www.teos-10.org
Other
23 stars 3 forks source link

gsw_f(lat) - (coriolis parameter) missing? #19

Closed willcoxe closed 1 year ago

willcoxe commented 2 years ago

I can't see the function referred to anywhere in the code, is it missing?

dankelley commented 2 years ago

Yes, it appears to be missing, since

git grep gsw_f\(

produces no output.

dankelley commented 2 years ago

I'm not the maintainer of GibbsSeaWater.jl, but perhaps that person will notice this and add it. In the meantime, you could define it in your own code, with

gsw_f(latitude) = 2*7.292115e-5*sin(latitude*pi/180.0)

(see https://www.teos-10.org/pubs/gsw/html/gsw_f.html)

willcoxe commented 2 years ago

Indeed. Just was surprised it was missing :) Thanks

Alexander-Barth commented 2 years ago

It is indeed surprising that this function was missing. As far as I know, the wrapper has been automatically generated by @kouketsu. But also the C library (that we are using) misses this function gsw_f:

$ nm /home/abarth/.julia/packages/GibbsSeaWater/vq3GE/deps/usr/lib/libgswteos-10.so | grep gsw_f
0000000000007650 T gsw_fdelta
000000000001a130 T gsw_frazil_properties
0000000000018b60 T gsw_frazil_properties_potential
000000000000dab0 T gsw_frazil_properties_potential_poly
000000000001a8f0 T gsw_frazil_ratios_adiabatic
000000000001b870 T gsw_frazil_ratios_adiabatic_poly
efiring commented 2 years ago

I imagine it was not included in the C library because it has nothing to do with seawater; it isn't used within gsw, and logically doesn't belong there.

kouketsu commented 2 years ago

As pointed out by @efiring, gsw_f is not included in GSW-C. I usually get the Coriolis parameter with the other package :).

I personally think GibbsSeaWater.jl doesn't need gsw_f (as @efiring pointing out).

dankelley commented 2 years ago

For context, GSW-R does not provide gsw_f(), because it works by accessing functions in the C source. Since the Julia version works similarly, I seems sensible to skip gsw_f(). I don't think there is much loss in this, because (a) the function is simple and has likely been part of user's general analysis code for a long time and (b) it does not relate to seawater.

Alexander-Barth commented 2 years ago

As this function is apparently in Matlab and python gsw package. I propose to include it here as well (as it is a single line) and then make a new release. Is this ok for you @kouketsu ? (I can make this change, but I am not in front of a computer right now). I would not be surprised if other people porting Matlab or python code to Julia and expecting this function to exist.

But if you prefer not to have it, it would be fine for me too.

kouketsu commented 2 years ago

I basically agree with @dankelley, but as pointed out @Alexander-Barth, I also understand that someone who plans to move Julia from Matlab may need gsw_f.

So I have added gsw_f, please try it (pkg> add GibbsSeaWater.jl#master; @Alexander-Barth and @willcoxe ).

willcoxe commented 2 years ago

I'd like to thank everyone who participated in the discussion. I never considered that the C code might not have implemented a function that is described in the teos-10 function contents (https://www.teos-10.org/pubs/gsw/html/gsw_contents.html) which is usually where I go to find my function descriptions irrespective of language (in my case python and julia).

I'm fine with it either way, but thanks for adding it!

Alexander-Barth commented 2 years ago

I tried (naively) latitude = 50 as an integer, and I got this error:

julia> GibbsSeaWater.gsw_f(50)
ERROR: InexactError: Int64(7.292115e-5)
Stacktrace:
 [1] Int64(x::Float64)
   @ Base ./float.jl:812
 [2] gsw_f(latitude::Int64)
   @ GibbsSeaWater ~/.julia/dev/GibbsSeaWater/src/utils.jl:1
 [3] top-level scope
   @ REPL[3]:1

I am wondering if we could use the following instead (julia's automatic type promotion):

gsw_f2(latitude) = 2 * 7.292115e-5 * sind(latitude) 
#(sind = sin for degrees)

But single precision floats become automatically double precision because we use a double precision constant. Maybe you wanted to avoid this?

What about this:

_Omega(::Type{Float32}) = 7.292115f-5  # constant as Float32
_Omega(::Type{T}) where T = 7.292115e-5 # default to Float64
gsw_f3(latitude::T) where T = 2 * _Omega(T) * sind(latitude)
julia> gsw_f3(50) # Int -> Float64
0.00011172168348669092

julia> gsw_f3(50.f0) # Float32 stays Float32
0.00011172169f0

julia> gsw_f3(50.0)
0.00011172168348669092

julia> gsw_f3(BigFloat(50.0)) # BigFloat works too :-)
0.0001117216834866909286569979820464075150286549620067559025634800693752850071411633

Sorry if this is overkill :-) I am likely to use only Float64, so the current master is fine for me (the C functions are Float64-only anyway as far as I know)

kouketsu commented 2 years ago

Thanks, @Alexander-Barth. It's my mistake, not overkill :). I fixed it.