Forums: » Register « |  Free Tools |  User CP |  Games |  Calendar |  Members |  FAQs |  Sitemap |  Support |
 User Name: Password: Remember me

New Free Tools on Dev Shed!
We're Excited to announce that Dev Shed now has 70 free tools on the site. To learn more, click here!

 Add This Thread To: Del.icio.us   Digg   Google   Spurl   Blink   Furl   Simpy   Y! MyWeb
 « Previous Thread | Next Thread » Thread Tools Search this Thread Rate Thread Display Modes
 Dev Shed Forums Sponsor:
#1
August 22nd, 2012, 03:34 PM
 newtoc2
Contributing User

Join Date: Aug 2012
Posts: 51
Time spent in forums: 15 h 33 m 10 sec
Reputation Power: 2
If/else decision-making for interpolation question

I've written a program that reads file data into 2 separate arrays.

Each line in the file contains 2 values: a and f(a)
So there are 2 arrays: first one is a, second one for f(a).
(In the question and code, array containing a is angle[N], and array containing f(a) is coeff[N].)

The user enters value b.

If b = a for any value in the first array, then no calculation is required, f(b) = f(a) and we print the element in the 2nd array containing with the same index as first array with b=a.

Else, we need 2 values from array containing a to interpolate f(b). One is larger than b, one smaller. And the corresponding f(a) and f(c) are used to get f(b)

Formula for interpolation:

f(b) = f(a) + ((b-a)/(c-a))*(f(c) - f(a));

I'm stuck on how to make it if b=a, then break out of the loop and just print f(a) and not print the f(b) computed with the formula. I tried using break and setting the counter k to the last count, but it still computes f(b) using the interpolation formula for all of the values accumulated before it had hit the point for b=a. (There's 2 values printed: one from the if statement where a=b, and one from the interpolation. It's confusing to the reader, the second one is also wrong.)

So I've resorted to the dumb "Ignore statement below" warning for the case b=a as to not confuse the user.

Any hints would be appreciated.

Code:
```/* Use linear interpolation to compute the coefficient of lift for flight-path angles. */
/* If user input is within bounds of file. */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define FILENAME "flight-path.txt"
#define PI 3.1415926535
#define N 17

int main(void)
{
/* Declare variables. */
double a, f_a, b, f_b, c, f_c;
double angle[N], coeff[N];
int k;
FILE *flight;

/* Open file. */
flight = fopen(FILENAME,"r");
if (flight == NULL)
printf("Error opening input file. \n");
else
{
/* Read data into two 1-D arrays. */
for (k=0;k<N;k++)
fscanf(flight,"%lf %lf",&angle[k],&coeff[k]);

/* Get user input from keyboard. */

printf("Enter new flight-path angle for interpolation \n"
"(between -4 and 21 degrees): \n");
scanf("%lf",&b);
printf(" b = %f \n",b);
a = -4;
c = 21;
f_a = -0.182;
f_c = 1.059;

/* Find f_a, a, f_c, c in arrays. */
for (k=0;k<=N-1;k++)
{
if (angle[k] == b)
{
printf("Data found in file: coefficient of lift at angle %f is %f. \n"
"Ignore statement below. \n",b,coeff[k]);
k = N-1;
}
if ((angle[k] <= b) && (angle[k] >= a))
{
a = angle[k];
f_a = coeff[k];
}

if ((angle[k] >= b) && (angle[k] <= c))
{
c = angle[k];
f_c = coeff[k];
}
}

/* Use linear interpolation to compute f_b. */
f_b = f_a + ((b-a)/(c-a))*(f_c - f_a);
b = b*PI/180;

/* Print new coefficient of lift. */
printf("New coefficient of lift at flight-path angle of %f radians: %f \n",b,f_b);
}
/* Exit Program. */
system("pause");
return 0;
}```

Thank you.

Last edited by newtoc2 : August 22nd, 2012 at 05:30 PM. Reason: If / else is a more accurate name

#2
August 22nd, 2012, 05:31 PM
 b49P23TIvg
Contributing User

Join Date: Aug 2011
Posts: 4,140
Time spent in forums: 1 Month 3 Weeks 2 Days 7 h 17 m 40 sec
Reputation Power: 455
1) why bother? The interpolant is valid at the ends of the intervals.

2) you could set an "already_computed" flag if you really must bypass the short recomputation.

3) learn about the "else" clause.

4) I'd separate the linear interpolation routine into a separate function to neaten your messy program. Linear interpolation is a good, reusable engineering tool.

Depending on how many times you will use your interpolant, it may pay off to precompute
(f_c - f_a) / (c - a)
for all the intervals.
__________________
[code]Code tags[/code] are essential for python code!

Last edited by b49P23TIvg : August 22nd, 2012 at 05:34 PM.

#3
August 23rd, 2012, 04:16 PM
 newtoc2
Contributing User

Join Date: Aug 2012
Posts: 51
Time spent in forums: 15 h 33 m 10 sec
Reputation Power: 2
1) why bother? The interpolant is valid at the ends of the intervals.

There's an error for division by 0. Since b=a=c, therefore c-a=0. The answer displayed is -1.#IND00, which I think says some irregularity occurred during the operation or some complex number showed up.

2) you could set an "already_computed" flag if you really must bypass the short recomputation.

3) learn about the "else" clause.

This is great!! New working code below. (I would've just posted the changed bits, but there's an else clause there so the braces are new. The code's not that long.)

Code:
```/* Use linear interpolation to compute the coefficient of lift for flight-path angles. */
/* If user input is within bounds of file. */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define FILENAME "flight-path.txt"
#define PI 3.1415926535
#define N 17

int main(void)
{
/* Declare variables. */
double a, f_a, b, f_b, c, f_c;
double angle[N], coeff[N], b_rad;
int k, already_computed;
FILE *flight;

/* Open file. */
flight = fopen(FILENAME,"r");
if (flight == NULL)
printf("Error opening input file. \n");
else
{
/* Read data into two 1-D arrays. */
for (k=0;k<N;k++)
fscanf(flight,"%lf %lf",&angle[k],&coeff[k]);

/* Get user input from keyboard. */

printf("Enter new flight-path angle for interpolation \n"
"(between -4 and 21 degrees): \n");
scanf("%lf",&b);
a = -4;
c = 21;
f_a = -0.182;
f_c = 1.059;

b_rad = b*PI/180;

/* Find f_a, a, f_c, c in arrays. */
for (k=0;k<=N-1;k++)
{
already_computed = (angle[k] == b);

if (already_computed)
{
printf("Found in data: coefficient of %f radians: %f \n",b_rad,coeff[k]);
break;
}
else
{
if ((angle[k] <= b) && (angle[k] >= a) && !already_computed)
{
a = angle[k];
f_a = coeff[k];
}

if ((angle[k] >= b) && (angle[k] <= c) && !already_computed)
{
c = angle[k];
f_c = coeff[k];
}
}
}

/* Use linear interpolation to compute f_b. */
f_b = f_a + ((b-a)/(c-a))*(f_c - f_a);

/* Print new coefficient of lift. */
printf("New coefficient of lift at flight-path angle of %f radians: %f \n",b_rad,f_b);
}
/* Exit Program. */
system("pause");
return 0;
}```

4) I'd separate the linear interpolation routine into a separate function to neaten your messy program. Linear interpolation is a good, reusable engineering tool.

Depending on how many times you will use your interpolant, it may pay off to precompute
(f_c - f_a) / (c - a)
for all the intervals.

Thanks for the tips!

#4
August 24th, 2012, 01:00 PM
 b49P23TIvg
Contributing User

Join Date: Aug 2011
Posts: 4,140
Time spent in forums: 1 Month 3 Weeks 2 Days 7 h 17 m 40 sec
Reputation Power: 455
math.h contains some definitions, for instance pi is M_PI . pi is wrong

If true in your code that a == c then either
1) the data is non-functional (multi-valued) and therefor invalid.
2) your interval search code is incorrect.

Interpolation (splines) routines usually validate their data. These routines require strictly increasing abscissa.

How do I handle non-functional data?
I parameterize it by integrated path length or by index number. Then I interpolate both x and y (and the other coordinates) with interpolants

x(t)
y(t)
etceteras

If in your code a and c are the same, do they have the same index? Does f(a) == f(c) ?

for instance:
Code:
```#define INRANGE(A,L,H) (((L)<=(A))&((A)<(H))) /* [ ) */

int strictly_increasing(float*a,size_t n) {
int i, result = 1;
for (i = 0; i < n-1; ++i)
result &= a[i] < a[i+1];
return result;
}

/* linearly_interpolate interpolates the discrete function given by f(x[0...n-1]) = y[0...n-1] at the point X , returning an approximation to y(X) */
float linearly_interpolate(float X,float*x,float*y,size_t n) {
int i;
char*func = "linearly_interpolate";/* c99  __func__ */
if (!n) {
fprintf(stderr,"\n%s: No data to interpolate, returning 0\n",func);
return 0;			/* or take stronger action */
}
if (1 == n)                   /* constant */
return y[0];
if (! strictly_increasing(x,n)) {
fprintf(stderr,"\n%s: Invalid abscissa, returning 0\n",func);
return 0;			/* or take stronger action */
}
for (i = 1; (i+1 < n) && (x[i] < X); ++i)
;
if (! INRANGE(X,x[i-1],x[i]))
fprintf(stderr,"\n%s:extrapolation\n",func); /* or take stronger action */
return y[i-1] + (X-x[i-1]) * ((y[i]-y[i-1])/(x[i]-x[i-1]));
}```

#5
August 27th, 2012, 06:57 AM
 newtoc2
Contributing User

Join Date: Aug 2012
Posts: 51
Time spent in forums: 15 h 33 m 10 sec
Reputation Power: 2
math.h contains some definitions, for instance pi is M_PI .

Thanks for that.

If true in your code that a == c then either
1) the data is non-functional (multi-valued) and therefor invalid.
2) your interval search code is incorrect.

Yes, my search code was incorrect -- I had [ ] instead of [ ) for the end limits.

Interpolation (splines) routines usually validate their data. These routines require strictly increasing abscissa.

Yes, my data had strictly increasing abscissa.

How do I handle non-functional data?
I parameterize it by integrated path length or by index number. Then I interpolate both x and y (and the other coordinates) with interpolants

Let me do a brief interp. of this code:

- INRANGE(a,b,c) means "a" is in interval [b,c)
- strictly_increasing is a function taking in an array and an integer specifying size of array. It returns a 0 or 1, indicating if the next value in the array is always greater than the previous.
- function linearly_interpolate takes in float X, which is the value we want to find f(X) of. It also takes in the array of x values, and the array of y values, plus an integer that specifies the size of the arrays
- if size is 0, then there's no values to compute. Return error msg.
- if size is 1, then there's only 1 value available, return that value.
- if the x coordinates are not strictly increasing, data is invalid. Return error message
- if X is not within bounds of the available x values in the array of x's, then we can extrapolate or just return error msg.
- else, we return the interpolated value. Our formulas are really the same, though syntax is different.

If in your code a and c are the same, do they have the same index? Does f(a) == f(c) ?

Yes. f(a) == f(c).

Thanks so much, you make some good points. I really appreciate it.

#6
August 27th, 2012, 07:38 AM
 b49P23TIvg
Contributing User

Join Date: Aug 2011
Posts: 4,140
Time spent in forums: 1 Month 3 Weeks 2 Days 7 h 17 m 40 sec
Reputation Power: 455
Our aren't quite the same: mine is reusable while you'll have to rewrite yours next time you need it.

#7
August 27th, 2012, 02:23 PM
 newtoc2
Contributing User

Join Date: Aug 2012
Posts: 51
Time spent in forums: 15 h 33 m 10 sec
Reputation Power: 2
By formulas are equal I meant

Formula 1:

f_b = f_a + ((b-a)/(c-a))*(f_c - f_a) ≡

f_b = f_a + ((b-a) * (f_c - f_a)) /(c-a) [By associativity of multiplication and division]

formula 2:

return y[i-1] + (X-x[i-1]) * ((y[i]-y[i-1])/(x[i]-x[i-1]))

syntax substitutes:

f_a = y[i-1]
b = X
a = x[i-1]
f_c = y[i]
c = x[i-1]

Code:
```/***** Not code, code tags used to preserve spacing. *****/

so y[i-1] + (X-x[i-1]) * ((y[i]-y[i-1])/(x[i]-x[i-1]))
= f_a    + (b - a)    * (f_c- f_a)    /(c - a)
≡ Fomula1```

I think you might be referring to how I set
a = -4;
c = 21;
f_a = -0.182;
f_c = 1.059;

so if I wanted to use a different set of data, then well -- have to reset those values or something more clever like your for loop. Of course, your code includes a much more comprehensive checking system as well. Kudos for that!

#8
August 27th, 2012, 03:17 PM
 b49P23TIvg
Contributing User

Join Date: Aug 2011
Posts: 4,140
Time spent in forums: 1 Month 3 Weeks 2 Days 7 h 17 m 40 sec
Reputation Power: 455
Actually, I was referring to that yours is in the function "main" with a project specific "coefficient of lift". If you need an interpolation routine with custom messages use the function interface and pass the message, possibly with varargs.

#9
August 29th, 2012, 01:30 PM
 newtoc2
Contributing User

Join Date: Aug 2012
Posts: 51
Time spent in forums: 15 h 33 m 10 sec
Reputation Power: 2
Quote:
 Originally Posted by SaraParker23 Best Luck for solving the issue..please let us know ,after solving issue

Actually, it's been solved already. Both the codes work -- the one with "already computed" predicate works but its redundant. Here's the version after I fixed the search limits from [] to [).

Code:
```		a = -4;
c = 21;
f_a = -0.182;
f_c = 1.059;

b_rad = b*PI/180;

/* Find f_a, a, f_c, c in arrays. */
for (k=0;k<=N-1;k++)
{
if ((angle[k] < b) && (angle[k] > a))
{
a = angle[k];
f_a = coeff[k];
}

if ((angle[k] > b) && (angle[k] < c))
{
c = angle[k];
f_c = coeff[k];
}
}
/* Use linear interpolation to compute f_b. */
f_b = f_a + ((b-a)/(c-a))*(f_c - f_a);

/* Print new coefficient of lift. */
printf("New coefficient of lift at flight-path angle of %f radians: %f \n",b_rad,f_b);
}
/* Exit Program. */
system("pause");
return 0;
}```

b49P23TIvg's code is more general and can be applied to more situations, also has more checking and error messages. (It's better is what I mean.)

I think we're just discussing what the differences are and what the pros and cons are for each code. But as you can see, his (assuming the Avatar is photo of author) method is making functions that can be referred to in any context, while mine is in the main function, no modularity, no reusability.
Comments on this post
b49P23TIvg agrees: Is me, few years ago. google earth near 44&deg; 8'55.36"N 73&deg;52'39.33"W Adirondack park, NY.
USA. with 24 hr. flu setting in.

 Viewing: Dev Shed Forums > Programming Languages > C Programming > If/then decision-making for interpolation question

## Developer Shed Advertisers and Affiliates

 Thread Tools Search this Thread Search this Thread: Advanced Search Display Modes Rate This Thread Linear Mode Rate This Thread: 5 : Excellent 4 : Good 3 : Average 2 : Bad 1 : Terrible

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts vB code is On Smilies are On [IMG] code is On HTML code is Off
 View Your Warnings | New Posts | Latest News | Latest Threads | Shoutbox Forum Jump Please select one User Control Panel Private Messages Subscriptions Who's Online Search Forums Forums Home -------------------- Programming Languages    PHP Development        PHP FAQs and Stickies    Perl Programming        Perl FAQs and Stickies    C Programming        C Programming FAQs and Stickies    Java Help        Java FAQs    Python Programming        Python Programming FAQs    Ruby Programming        Ruby Programming FAQs    Game Development        Game Development FAQs Programming Languages - More    ASP Programming        ASP Programming FAQs    .Net Development        .Net Development FAQs    Visual Basic Programming        Visual Basic Programming FAQs    Software Design        Software Design FAQs    ColdFusion Development        ColdFusion Development FAQs    Delphi Programming        Delphi Programming FAQs    Regex Programming        Regex Programming FAQs    XML Programming        XML Programming FAQs    Other Programming Languages        Other Programming Languages FAQs Web Design    HTML Programming        HTML Programming FAQs    JavaScript Development        JavaScript Development FAQs    CSS Help        CSS Help FAQs    Flash Help        Flash Help FAQs    Photoshop Help        Photoshop Help FAQs    Web Design Help        Web Design Help FAQs    Website Critiques        Website Critiques FAQs    Search Engine Optimization        Search Engine Optimization FAQs Mobile Programming    Mobile Programming        Mobile Programming FAQs    iPhone SDK Development        iPhone SDK Development FAQs    Android Development        Android Development FAQs    BlackBerry Development        BlackBerry Development FAQs Web Site Management    Business Help        Business Help FAQs    Development Software        Development Software FAQs    Scripts        Scripts FAQs Databases    Database Management        Database Management FAQs    DB2 Development        DB2 Development FAQs    MySQL Help        MySQL Help FAQs    PostgreSQL Help        PostgreSQL Help FAQs    Firebird SQL Development        Firebird SQL Development FAQs    MS SQL Development        MS SQL Development FAQs    Oracle Development        Oracle Development FAQs    LDAP Programming        LDAP Programming FAQs System Administration    Mail Server Help        Mail Server Help FAQs    Apache Development        Apache Development FAQs    Security and Cryptography        Security and Cryptography FAQs    Antivirus Protection        Antivirus Protection FAQs    DNS        DNS FAQs    IIS        IIS FAQs    Networking Help        Networking Help FAQs    FTP Help        FTP Help FAQs Operating Systems    BSD Help        BSD Help FAQs    Linux Help        Linux Help FAQs    UNIX Help        UNIX Help FAQs    Windows Help        Windows Help FAQs    Mac Help        Mac Help FAQs Web Hosting    Web Hosting        Web Hosting FAQs    Free Web Hosting        Free Web Hosting FAQs    Web Hosting Requests        Web Hosting Requests FAQs    Web Hosting Offers        Web Hosting Offers FAQs Computer Hardware    Computer Hardware    CPUs        CPUs FAQs    Cooling        Cooling FAQs    Embedded Programming        Embedded Programming FAQs    Motherboards        Motherboards FAQs    Multimedia Hardware        Multimedia Hardware FAQs Other    Dev Shed Lounge        Dev Shed Lounge FAQs    Development Articles        Development Articles FAQs    Beginner Programming        Beginner Programming FAQs    Hire A Programmer        Hire A Programmer FAQs    Project Help Wanted        Project Help Wanted FAQs Latest News Updated Hourly    Technology News    Business News    Science News Forum Information    Forum Rules/Guidelines        Forum Rules/Guidelines FAQs    Forum Announcements        Forum Announcements FAQs    Dev Shed Gaming Center        Go to the Dev Shed Battle Arena        Go to the Dev Shed Arcade Games        Go to the Legend of the Green Dragon    Suggestions & Feedback        Suggestions & Feedback FAQs

 Forums: » Register « |  Free Tools |  User CP |  Games |  Calendar |  Members |  FAQs |  Sitemap |  Support |

Powered by: vBulletin Version 3.0.5
Copyright ©2000 - 2013, Jelsoft Enterprises Ltd.

© 2003-2013 by Developer Shed. All rights reserved. DS Cluster - Follow our Sitemap