source: issm/oecreview/Archive/13393-13976/ISSM-13568-13569.diff@ 13980

Last change on this file since 13980 was 13980, checked in by Mathieu Morlighem, 12 years ago

preparing oecreview for 13393-13976'
File size: 5.8 KB
RevLine 
[13980]1Index: ../trunk-jpl/src/c/shared/Sorting/binary_search.cpp
2===================================================================
3--- ../trunk-jpl/src/c/shared/Sorting/binary_search.cpp (revision 13568)
4+++ ../trunk-jpl/src/c/shared/Sorting/binary_search.cpp (revision 13569)
5@@ -10,7 +10,7 @@
6
7 #include <stdio.h>
8
9-int binary_search(int* poffset,int target,int* sorted_integers,int num_integers){
10+int binary_search(int* poffset,int target,int* sorted_integers,int num_integers){ /*{{{*/
11
12 /*output: */
13 int offset; //offset, if found
14@@ -62,4 +62,68 @@
15
16 /*Return result: */
17 return found;
18-}
19+} /*}}}*/
20+int binary_search(int* poffset,double target,double* sorted_doubles,int num_doubles){ /*{{{*/
21+
22+ /*output: */
23+ int offset=0;
24+ int found=0; /*found=0: not found.
25+ found=1: found, and target is == to value at offset
26+ found=2: found, and target is > to value at offset and < to value at offset+1
27+ */
28+
29+ /*intermediary: */
30+ double* beg=NULL;
31+ double* end=NULL;
32+ double* mid=NULL;
33+
34+ // point to beginning and end of the array
35+ beg=sorted_doubles;
36+ end=sorted_doubles+num_doubles;
37+ mid=beg+(int)(num_doubles/2.0);
38+
39+ if (target<*beg){
40+ offset==-1;
41+ found=0;
42+ }
43+ if (*beg==target){
44+ found=1;
45+ offset=0;
46+ }
47+ else if(*(end-1)==target){
48+ found=1;
49+ offset=num_doubles-1;
50+ }
51+ else{
52+ while((beg <= end) && !( target>=*mid && target<*(mid+1)) ){
53+ // is the target in lower or upper half?
54+ if (target < *mid) {
55+ end = mid - 1; //new end
56+ mid = beg + (end-beg)/2; //new middle
57+ }
58+ else {
59+ beg = mid + 1; //new beginning
60+ mid = beg + (end-beg)/2; //new middle
61+ }
62+ }
63+
64+ //did we find the target?
65+ if( target>*mid && target<*(mid+1)){
66+ found=2;
67+ offset=mid-sorted_doubles;
68+ }
69+ else if( target==*mid){
70+ found=1;
71+ offset=mid-sorted_doubles;
72+ }
73+ else {
74+ found=0;
75+ }
76+ }
77+
78+ /*Assign output pointers:*/
79+ *poffset=offset;
80+
81+ /*Return result: */
82+ return found;
83+} /*}}}*/
84Index: ../trunk-jpl/src/c/shared/Sorting/sorting.h
85===================================================================
86--- ../trunk-jpl/src/c/shared/Sorting/sorting.h (revision 13568)
87+++ ../trunk-jpl/src/c/shared/Sorting/sorting.h (revision 13569)
88@@ -6,5 +6,6 @@
89 #define _SORTING_H_
90
91 int binary_search(int* poffset,int target,int* sorted_integers,int num_integers);
92+int binary_search(int* poffset,double target,double* sorted_doubles,int num_doubles);
93
94 #endif //ifndef _SORTING_H_
95Index: ../trunk-jpl/src/c/classes/objects/Inputs/TransientInput.cpp
96===================================================================
97--- ../trunk-jpl/src/c/classes/objects/Inputs/TransientInput.cpp (revision 13568)
98+++ ../trunk-jpl/src/c/classes/objects/Inputs/TransientInput.cpp (revision 13569)
99@@ -425,56 +425,58 @@
100 /*FUNCTION TransientInput::GetTimeInput{{{*/
101 Input* TransientInput::GetTimeInput(IssmDouble intime){
102
103- int i,j;
104 IssmDouble deltat;
105 IssmDouble alpha1,alpha2;
106- bool found=false;
107+
108 Input* input=NULL;
109 Input* input1=NULL;
110 Input* input2=NULL;
111+
112+ bool found=false;
113+ int found_offset=0;
114+ int offset;
115
116 /*Ok, we have the time, go through the timesteps, and figure out which interval we
117 *fall within. Then interpolate the values on this interval: */
118- if(intime<this->timesteps[0]){
119- /*get values for the first time: */
120- input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy();
121- found=true;
122+ found_offset=binary_search(&offset,intime,this->timesteps,this->numtimesteps);
123+
124+ if (found_offset==0){
125+ if(intime<this->timesteps[0]){
126+ /*get values for the first time: */
127+ input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy();
128+ found=true;
129+ }
130+ else if(intime>this->timesteps[this->numtimesteps-1]){
131+ /*get values for the last time: */
132+ input=(Input*)((Input*)this->inputs->GetObjectByOffset(numtimesteps-1))->copy();
133+ found=true;
134+ }
135 }
136- else if(intime>this->timesteps[this->numtimesteps-1]){
137- /*get values for the last time: */
138- input=(Input*)((Input*)this->inputs->GetObjectByOffset(numtimesteps-1))->copy();
139- found=true;
140- }
141 else{
142- /*Find which interval we fall within: */
143- for(i=0;i<this->numtimesteps;i++){
144- if(intime==this->timesteps[i]){
145- /*We are right on one step time: */
146- input=(Input*)((Input*)this->inputs->GetObjectByOffset(i))->copy();
147- found=true;
148- break; //we are done with the time interpolation.
149- }
150- else{
151- if(this->timesteps[i]<intime && intime<this->timesteps[i+1]){
152- /*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */
153- deltat=this->timesteps[i+1]-this->timesteps[i];
154- alpha2=(intime-this->timesteps[i])/deltat;
155- alpha1=(1.0-alpha2);
156+ if (found_offset==1){
157+ /*intime is exactly equal to this->timesteps[offset]: */
158+ input=(Input*)((Input*)this->inputs->GetObjectByOffset(offset))->copy();
159+ found=true;
160+ }
161+ else if (found_offset==2){
162+
163+ /*ok, we have the interval ]offset:offset+1[. Interpolate linearly for now: */
164+ deltat=this->timesteps[offset+1]-this->timesteps[offset];
165+ alpha2=(intime-this->timesteps[offset])/deltat;
166+ alpha1=(1.0-alpha2);
167
168- input1=(Input*)this->inputs->GetObjectByOffset(i);
169- input2=(Input*)this->inputs->GetObjectByOffset(i+1);
170+ input1=(Input*)this->inputs->GetObjectByOffset(offset);
171+ input2=(Input*)this->inputs->GetObjectByOffset(offset+1);
172
173- input=(Input*)input1->copy();
174- input->Scale(alpha1);
175- input->AXPY(input2,alpha2);
176+ input=(Input*)input1->copy();
177+ input->Scale(alpha1);
178+ input->AXPY(input2,alpha2);
179
180- found=true;
181- break;
182- }
183- else continue; //keep looking on the next interval
184- }
185+ found=true;
186 }
187+ else _error_("binary_search returned the following value for found: " << found_offset<< "which is not supported yet!");
188 }
189+
190 if(!found)_error_("did not find time interval on which to interpolate forcing values!");
191
192 /*Assign output pointer*/
Note: See TracBrowser for help on using the repository browser.