#1
  1. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Aug 2012
    Posts
    51
    Rep Power
    3

    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 04:30 PM. Reason: If / else is a more accurate name
  2. #2
  3. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,895
    Rep Power
    481
    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.
    Last edited by b49P23TIvg; August 22nd, 2012 at 04:34 PM.
    [code]Code tags[/code] are essential for python code and Makefiles!
  4. #3
  5. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Aug 2012
    Posts
    51
    Rep Power
    3
    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!
  6. #4
  7. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,895
    Rep Power
    481
    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]));
    }
    [code]Code tags[/code] are essential for python code and Makefiles!
  8. #5
  9. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Aug 2012
    Posts
    51
    Rep Power
    3
    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.
  10. #6
  11. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,895
    Rep Power
    481
    Our aren't quite the same: mine is reusable while you'll have to rewrite yours next time you need it.
    [code]Code tags[/code] are essential for python code and Makefiles!
  12. #7
  13. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Aug 2012
    Posts
    51
    Rep Power
    3
    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!
  14. #8
  15. Contributing User
    Devshed Demi-God (4500 - 4999 posts)

    Join Date
    Aug 2011
    Posts
    4,895
    Rep Power
    481
    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.
    [code]Code tags[/code] are essential for python code and Makefiles!
  16. #9
  17. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Aug 2012
    Posts
    51
    Rep Power
    3
    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 8'55.36"N 7352'39.33"W Adirondack park, NY. USA. with 24 hr. flu setting in.

IMN logo majestic logo threadwatch logo seochat tools logo