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 }