001/**
002 *
003 * Copyright (c) 2014, the Railo Company Ltd. All rights reserved.
004 *
005 * This library is free software; you can redistribute it and/or
006 * modify it under the terms of the GNU Lesser General Public
007 * License as published by the Free Software Foundation; either 
008 * version 2.1 of the License, or (at your option) any later version.
009 * 
010 * This library is distributed in the hope that it will be useful,
011 * but WITHOUT ANY WARRANTY; without even the implied warranty of
012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
013 * Lesser General Public License for more details.
014 * 
015 * You should have received a copy of the GNU Lesser General Public 
016 * License along with this library.  If not, see <http://www.gnu.org/licenses/>.
017 * 
018 **/
019package lucee.runtime.img.interpolation;
020
021public class Bessel implements Interpolation
022{
023    private double J1(double x) {
024        double[] Pone
025            = { 5.811993540016061E20, -6.672106568924916E19,
026                2.3164335806340024E18, -3.588817569910106E16,
027                2.9087952638347756E14, -1.3229834803321265E12,
028                3.4132341823017006E9, -4695753.530642996, 2701.1227108923235 };
029        double[] Qone = { 1.1623987080032122E21, 1.185770712190321E19,
030                          6.0920613989175216E16, 2.0816612213076075E14,
031                          5.2437102621676495E11, 1.013863514358674E9,
032                          1501793.5949985855, 1606.9315734814877, 1.0 };
033        double p = Pone[8];
034        double q = Qone[8];
035        for (int i = 7; i >= 0; i--) {
036            p = p * x * x + Pone[i];
037            q = q * x * x + Qone[i];
038        }
039        return p / q;
040    }
041    
042    private double P1(double x) {
043        double[] Pone
044            = { 35224.66491336798, 62758.84524716128, 31353.963110915956,
045                4985.4832060594335, 211.15291828539623, 1.2571716929145342 };
046        double[] Qone
047            = { 35224.66491336798, 62694.34695935605, 31240.406381904104,
048                4930.396490181089, 203.07751891347593, 1.0 };
049        double p = Pone[5];
050        double q = Qone[5];
051        for (int i = 4; i >= 0; i--) {
052            p = p * (8.0 / x) * (8.0 / x) + Pone[i];
053            q = q * (8.0 / x) * (8.0 / x) + Qone[i];
054        }
055        return p / q;
056    }
057    
058    private double Q1(double x) {
059        double[] Pone
060            = { 351.17519143035526, 721.0391804904475, 425.98730116544425,
061                83.18989576738508, 4.568171629551227, 0.03532840052740124 };
062        double[] Qone
063            = { 7491.737417180912, 15414.177339265098, 9152.231701516992,
064                1811.1867005523513, 103.81875854621337, 1.0 };
065        double p = Pone[5];
066        double q = Qone[5];
067        for (int i = 4; i >= 0; i--) {
068            p = p * (8.0 / x) * (8.0 / x) + Pone[i];
069            q = q * (8.0 / x) * (8.0 / x) + Qone[i];
070        }
071        return p / q;
072    }
073    
074    private double BesselOrderOne(double x) {
075        if (x == 0.0)
076            return 0.0;
077        double p = x;
078        if (x < 0.0)
079            x = -x;
080        if (x < 8.0)
081            return p * J1(x);
082        double q
083            = (Math.sqrt(2.0 / (3.141592653589793 * x))
084               * (P1(x) * (1.0 / Math.sqrt(2.0) * (Math.sin(x) - Math.cos(x)))
085                  - 8.0 / x * Q1(x) * (-1.0 / Math.sqrt(2.0)
086                                       * (Math.sin(x) + Math.cos(x)))));
087        if (p < 0.0)
088            q = -q;
089        return q;
090    }
091    
092    public double f(double x) {
093        if (x == 0.0)
094            return 0.7853981633974483;
095        return BesselOrderOne(3.141592653589793 * x) / (2.0 * x);
096    }
097    
098    public double getSupport() {
099        return 3.2383;
100    }
101}