source: trunk/src/HuntLocator.cpp@ 2729

Last change on this file since 2729 was 2727, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: Yes CAS-4770, CAS-4774

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Updated STApplyCal to be able to specify interpolation method.
The method can be specified in time and frequency axes independently.
Possible options are nearest, linear (default), (natural) cubic spline,
and polynomial with arbitrary order.

File size: 1.8 KB
RevLine 
[2727]1//
2// C++ Implementation: HuntLocator
3//
4// Description:
5//
6//
7// Author: Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2012
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <assert.h>
13
14#include "HuntLocator.h"
15
16namespace asap {
17
18HuntLocator::HuntLocator(double *v, unsigned int n)
19 : Locator(v, n),
20 prev_(0)
21{}
22
23HuntLocator::~HuntLocator()
24{}
25
26unsigned int HuntLocator::locate(double x)
27{
28 if (n_ == 1)
29 return 0;
30 bool ascending = (x_[n_-1] >= x_[0]);
31 if (ascending) {
32 if (x <= x_[0])
33 return 0;
34 else if (x > x_[n_-1])
35 return n_;
36 }
37 else {
38 if (x > x_[0])
39 return 0;
40 else if (x <= x_[n_-1])
41 return n_;
42 }
43
44 unsigned int jl = 0;
45 unsigned int ju = n_;
46 if (prev_ > 0 && prev_ < n_) {
47 // do hunt
48 jl = prev_;
49 unsigned int inc = 1;
50 if ((x >= x_[jl]) == ascending) {
51 // hunt up
52 if (jl >= n_ - 1)
53 return jl;
54 ju = jl + inc;
55 while ((x >= x_[ju]) == ascending) {
56 jl = ju;
57 inc <<= 1;
58 ju = jl + inc;
59 if (ju > n_ - 1) {
60 ju = n_;
61 break;
62 }
63 }
64 }
65 else {
66 // hunt down
67 if (jl == 0)
68 return jl;
69 ju = jl;
70 jl -= inc;
71 while ((x < x_[jl]) == ascending) {
72 ju = jl;
73 inc <<= 1;
74 if (inc >= ju) {
75 jl = 0;
76 break;
77 }
78 else
79 jl = ju - inc;
80 }
81 }
82 }
83
84 // final bisection phase
85 unsigned int jm;
86 if (ascending) {
87 while (ju - jl > 1) {
88 jm = (ju + jl) >> 1;
89 if (x > x_[jm])
90 jl = jm;
91 else
92 ju = jm;
93 }
94 }
95 else {
96 while (ju - jl > 1) {
97 jm = (ju + jl) >> 1;
98 if (x < x_[jm])
99 jl = jm;
100 else
101 ju = jm;
102 }
103 }
104 prev_ = jl;
105 return ju;
106}
107
108}
Note: See TracBrowser for help on using the repository browser.