001    package railo.runtime.img.interpolation;
002    
003    public class Bessel implements Interpolation
004    {
005        private double J1(double x) {
006            double[] Pone
007                = { 5.811993540016061E20, -6.672106568924916E19,
008                    2.3164335806340024E18, -3.588817569910106E16,
009                    2.9087952638347756E14, -1.3229834803321265E12,
010                    3.4132341823017006E9, -4695753.530642996, 2701.1227108923235 };
011            double[] Qone = { 1.1623987080032122E21, 1.185770712190321E19,
012                              6.0920613989175216E16, 2.0816612213076075E14,
013                              5.2437102621676495E11, 1.013863514358674E9,
014                              1501793.5949985855, 1606.9315734814877, 1.0 };
015            double p = Pone[8];
016            double q = Qone[8];
017            for (int i = 7; i >= 0; i--) {
018                p = p * x * x + Pone[i];
019                q = q * x * x + Qone[i];
020            }
021            return p / q;
022        }
023        
024        private double P1(double x) {
025            double[] Pone
026                = { 35224.66491336798, 62758.84524716128, 31353.963110915956,
027                    4985.4832060594335, 211.15291828539623, 1.2571716929145342 };
028            double[] Qone
029                = { 35224.66491336798, 62694.34695935605, 31240.406381904104,
030                    4930.396490181089, 203.07751891347593, 1.0 };
031            double p = Pone[5];
032            double q = Qone[5];
033            for (int i = 4; i >= 0; i--) {
034                p = p * (8.0 / x) * (8.0 / x) + Pone[i];
035                q = q * (8.0 / x) * (8.0 / x) + Qone[i];
036            }
037            return p / q;
038        }
039        
040        private double Q1(double x) {
041            double[] Pone
042                = { 351.17519143035526, 721.0391804904475, 425.98730116544425,
043                    83.18989576738508, 4.568171629551227, 0.03532840052740124 };
044            double[] Qone
045                = { 7491.737417180912, 15414.177339265098, 9152.231701516992,
046                    1811.1867005523513, 103.81875854621337, 1.0 };
047            double p = Pone[5];
048            double q = Qone[5];
049            for (int i = 4; i >= 0; i--) {
050                p = p * (8.0 / x) * (8.0 / x) + Pone[i];
051                q = q * (8.0 / x) * (8.0 / x) + Qone[i];
052            }
053            return p / q;
054        }
055        
056        private double BesselOrderOne(double x) {
057            if (x == 0.0)
058                return 0.0;
059            double p = x;
060            if (x < 0.0)
061                x = -x;
062            if (x < 8.0)
063                return p * J1(x);
064            double q
065                = (Math.sqrt(2.0 / (3.141592653589793 * x))
066                   * (P1(x) * (1.0 / Math.sqrt(2.0) * (Math.sin(x) - Math.cos(x)))
067                      - 8.0 / x * Q1(x) * (-1.0 / Math.sqrt(2.0)
068                                           * (Math.sin(x) + Math.cos(x)))));
069            if (p < 0.0)
070                q = -q;
071            return q;
072        }
073        
074        public double f(double x) {
075            if (x == 0.0)
076                return 0.7853981633974483;
077            return BesselOrderOne(3.141592653589793 * x) / (2.0 * x);
078        }
079        
080        public double getSupport() {
081            return 3.2383;
082        }
083    }