Refactored obsolescent code

[jacobwilliams]

Fixes #19

[certik]

Thanks for doing that. I left some minor comments, that perhaps make the code even simpler / more readable.

My main comment is this: does the new code actually work? It's very easy to make a tiny mistake in such refactoring (even though I haven't found any by just looking at it), which will make the code fail in various ways.

One reason I just took the subroutines verbatim from scipy is that they were extensively tested in scipy.

Given that the changes are essentially just syntax changes, it does seem it should work. So maybe it's fine.

Did you test the subroutines to ensure they work as before? If so, then I think we can merge it.

[certik]

You can also write this as 1._dp / 3.

[certik]

You can also just assign 0, it will be converted to double precision.

[certik]

One can just do sj(0) = 1.

[jacobwilliams]

Agreed. I'm usually paranoid about the kinds, so I usually keep the _dp even for 0 and 1.

[certik]

Since f1 is double precision, this can also be written as:

f=(2*k+3)*f1/x-f0

[certik]

The original is F1=1.0D0-100, which is a bit confusing, but I think it's just 1-100 which is -99. Since f1 is double precision, you can just write it as 1-100.

[jacobwilliams]

Yes, I noticed that. I actually wondered if that was some kind of typo (is it possible they meant 1.0e-100?)

[certik]

I first thought you made a typo, but then the original is the same. I wonder that too, it really looks like a typo. I don't remember with these, but I know it's possible that some of these old Fortran routines can produce incorrect answers, a different example that I found is:

https://mail.scipy.org/pipermail/scipy-user/2012-September/033189.html

(This has since been fixed.)

[certik]

dj(k)=sj(k-1)-(k+1)*sj(k)/x

[jacobwilliams]

Test case (same result from refactored version and the one in scipy):

 n    =           2
 mp   =           5
 x    =   1.2345600000000001     
 sj   =  0.76464742242968342       0.35211808607480077        9.1005002907604604E-002
 dj   = -0.35211808607480077       0.19421247220482474       0.13097451369045834     
 nm1  =           2
 sy   = -0.26725029005084994      -0.98112154280524222       -2.1168903174657765     
 dy   =  0.98112154280524222        1.3221767816269010        4.1629548507255132     
 nm2  =           2
 msta1=           6
 msta2=          16

This is the only one I did.

[jacobwilliams]

Here's another one (same result for both):

 n    =           4
 mp   =           6
 x    =   10.234560000000000     
 sj   =  -7.0754079663227554E-002   6.0471793911910328E-002   8.8479842347283441E-002  -1.7245782071985880E-002 -0.10027521747859354     
 dj   =  -6.0471793911910390E-002  -8.2571254785931469E-002   3.4536186807955657E-002   9.5220056708032075E-002   3.1742751621008022E-002
 nm1  =           4
 sy   =   6.7385044668486882E-002   7.7338148218054237E-002  -4.4715340777536736E-002  -9.9183416005597977E-002  -2.3121861020995477E-002
 dy   =  -7.7338148218054237E-002   5.2271908741186784E-002   9.0445308890580486E-002  -5.9512254640897572E-003  -8.7887443818715744E-002
 nm2  =           4
 msta1=          22
 msta2=          32

[certik]

@jacobwilliams thanks for doing the tests. If you wouldn't mind doing a few more just to make sure, that would be great. If there is a bug introduced by this, we can always trace it to this PR and fix it.