Some (probably) ridiculous questions

Some (probably) ridiculous questions

[Jean-Max Redonnet]

Sorry, if this point is obvious, but I'm not very familiar with GSL and I 
would like to use it in the right way.

My question is : How to declare a function of multiple variables ?

I manage with parametric curves and surfaces, so I need functions of one or 
two parameters. Futhermore This functions should be vectorial.
For example : I need to deal with a ruled surface.
Its definition is
C0(u) : vectorial function with 3 components : C0x(u), C0y(u) and C0z(u)
C1(u) : idem with C1x(u), C1y(u) and C1z(u)
then surface S(u,v) is built from C0 and C1 with
S(u,v) = (1-v)C0(u) + vC1(u)

How can I define this with GSL ?

For moment I've defined a
  struct v3func
    {
      gsl_function fx;
      gsl_function fy;
      gsl_function fz;
    };

then, once declared
  struct v3func *C0;
I can do
  C0 = (struct v3func *)malloc(sizeof(struct v3func *));
  C0->fx.function = &C0_fx;
  C0->fy.function = &C0_fy;
  C0->fz.function = &C0_fz;
with C0_fx, C0_fy and C0_fz declared by
  double C0_fx(double u, void * params);
  double C0_fy(double u, void * params);
  double C0_fz(double u, void * params);

This seems to work. Idem for C1 (of course).

But I have a problem to declare S(u,v)...
I've tried to do the same with S_fx(double u, double v, void * params)
but compilation says : "assignment from incompatible pointer type"

Any help would be very appreciated...

-- 
===============================================================================
Jean-Max Redonnet
PhD University Paul Sabatier - Toulouse (France)
	mailto:redonnetNO@SPAMlgmt.ups-tlse.fr

Re: Some (probably) ridiculous questions

[Brian Gough]

Jean-Max Redonnet writes:
 > Sorry, if this point is obvious, but I'm not very familiar with GSL and I 
 > would like to use it in the right way.
 > 
 > My question is : How to declare a function of multiple variables ?
 > 
 > I manage with parametric curves and surfaces, so I need functions of one or 
 > two parameters. Futhermore This functions should be vectorial.
 > For example : I need to deal with a ruled surface.

Hello,

The gsl_function type is mainly designed for interfacing to the
univariate gsl routines.  It is not really a general facility so it
does not handle other cases that are not used in the library, such as
surfaces.

The way to handle this sort of situation is,

-- define your own surface function type in the way that you normally
would in your C-programs

-- if you need to use a gsl routine, for example to integrate along a
line, define a function to create the appropriate gsl_function by
mapping from your surface type.

Philosophically, calls to GSL routines, and the GSL types, can be at a
"lower-level" than your program.

regards
Brian Gough

Re: Some (probably) ridiculous questions

[Jean-Max Redonnet]

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1: Type: text/plain, Size: 2690 bytes --]

Le Vendredi  6 Juillet 2001 17:35, Brian Gough a écrit :

<snip>

>
> Hello,
>
> The gsl_function type is mainly designed for interfacing to the
> univariate gsl routines.  It is not really a general facility so it
> does not handle other cases that are not used in the library, such as
> surfaces.
>
> The way to handle this sort of situation is,
>
> -- define your own surface function type in the way that you normally
> would in your C-programs
>
> -- if you need to use a gsl routine, for example to integrate along a
> line, define a function to create the appropriate gsl_function by
> mapping from your surface type.
>
> Philosophically, calls to GSL routines, and the GSL types, can be at a
> "lower-level" than your program.
>
> regards
> Brian Gough

Thanks for your help to a definitively poor programmer (and mathematician)

But I'm not sure I really understand what you mean...

I think I have defined a surface function type as you proposed.
This part seems to work as expected.
but I don't see how to define a function to create the appropriate 
gsl_function...

Here the code :

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_diff.h>

#define RAD2DEG(X) (X*180/M_PI)
#define DEG2RAD(X) (X*M_PI/180)

struct surf_function_struct {
  double (* fx)(double u, double v);
  double (* fy)(double u, double v);
  double (* fz)(double u, double v);;
};

typedef struct surf_function_struct surf_function ;

double C0_fx(double u)
{
  return 80*u;
}

double C0_fy(double u)
{
  return 80*tan(DEG2RAD(22.5))*(u-0.5);
}

double C0_fz(double u)
{
  return 40;
}

double C1_fx(double u)
{
  return 80*u;
}

double C1_fy(double u)
{
  return -80*tan(DEG2RAD(22.5))*(u-0.5);
}

double C1_fz(double u)
{
  return -40;
}

double S_fx(double u, double v)
{
  return (1-v)*C0_fx(u)+v*C1_fx(u);
}

double S_fy(double u, double v)
{
  return (1-v)*C0_fy(u)+v*C1_fy(u);
}

double S_fz(double u, double v)
{
  return (1-v)*C0_fz(u)+v*C1_fz(u);
}

int main (void)
{
  double u,v;
  surf_function surf;

  u = 0.5;
  v = 0.5;
  surf.fx = &S_fx;
  surf.fy = &S_fy;
  surf.fz = &S_fz;
  
  printf("              | %f \n", S_fx(u,v));
  printf(" S(%.1f,%.1f) = | %f \n", u, v, S_fy(u,v));
  printf("              | %f \n", S_fz(u,v));

  printf("              | %f \n", surf.fx(u,v));
  printf(" S(%.1f,%.1f) = | %f \n", u, v, surf.fy(u,v));
  printf("              | %f \n", surf.fz(u,v));
  
  return 0;
}

You would be my definitive heroe if you can help me on this point.

Thanks.
-- 
===============================================================================
Veuillez noter ma nouvelle adresse électronique :
	mailto:redonnetNO@SPAMlgmt.ups-tlse.fr

Re: Some (probably) ridiculous questions

[Brian Gough]

Jean-Max Redonnet writes:
 > I think I have defined a surface function type as you proposed.
 > This part seems to work as expected.
 > but I don't see how to define a function to create the appropriate 
 > gsl_function...
 > 

The gsl_function has to be a parametric function.  The example below
defines a gsl_function 'line.f' for |F|^2 on a path P(t),

  (u(t),v(t)) = (0.5,0.5) + 0.1*(cos(t),sin(t))

and then finds the minimum (which is at t=pi) using the standard gsl
minimization routines.  The relevant functions are line_eval() and
make_line_function().

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_diff.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_errno.h>

#define RAD2DEG(X) (X*180/M_PI)
#define DEG2RAD(X) (X*M_PI/180)

struct surf_function_struct {
  double (* fx)(double u, double v);
  double (* fy)(double u, double v);
  double (* fz)(double u, double v);;
};

typedef struct surf_function_struct surf_function ;

double C0_fx(double u)
{
  return 80*u;
}

double C0_fy(double u)
{
  return 80*tan(DEG2RAD(22.5))*(u-0.5);
}

double C0_fz(double u)
{
  return 40;
}

double C1_fx(double u)
{
  return 80*u;
}

double C1_fy(double u)
{
  return -80*tan(DEG2RAD(22.5))*(u-0.5);
}

double C1_fz(double u)
{
  return -40;
}

double S_fx(double u, double v)
{
  return (1-v)*C0_fx(u)+v*C1_fx(u);
}

double S_fy(double u, double v)
{
  return (1-v)*C0_fy(u)+v*C1_fy(u);
}

double S_fz(double u, double v)
{
  return (1-v)*C0_fz(u)+v*C1_fz(u);
}


struct path_function_struct {
  void (* uv)(double t, double *u, double *v);
};

typedef struct path_function_struct path_function ;

void P_uv(double t, double *u, double *v)
{
  *u = 0.5 + 0.1 * cos(t);
  *v = 0.5 + 0.1 * sin(t);
}

struct line_function_struct
{
  path_function p;
  surf_function s;
  gsl_function f;
};

typedef struct line_function_struct line_function ;

double 
line_eval (double t, void * params)
{
  line_function * line = (line_function *) params;

  double u, v;

  line->p.uv(t, &u, &v);
  
  {
    double a = line->s.fx(u,v);
    double b = line->s.fy(u,v);
    double c = line->s.fz(u,v);

    /* compute |F|^2 */

    return a*a + b*b + c*c;
  }
}

void
make_line_function (path_function path, surf_function surf, 
		    line_function * line)
{
  line->p = path;
  line->s = surf;
  line->f.function = &line_eval;
  line->f.params = line;
}


int main (void)
{
  double u,v,t;
  surf_function surf;
  path_function path;

  surf.fx = &S_fx;
  surf.fy = &S_fy;
  surf.fz = &S_fz;
  
  path.uv = &P_uv;

  for (t = 0 ; t < 2*M_PI; t+= 0.5)
    {
      path.uv(t, &u, &v);

      printf("                  | %f \n", S_fx(u,v));
      printf(" S(%.3f,%.3f) = | %f \n", u, v, S_fy(u,v));
      printf("                  | %f \n", S_fz(u,v));
      
      printf("                  | %f \n", surf.fx(u,v));
      printf(" S(%.3f,%.3f) = | %f \n", u, v, surf.fy(u,v));
      printf("                  | %f \n", surf.fz(u,v));
    }

  {
    line_function line;

    make_line_function(path, surf, &line);

    {
      int status, iterations = 0;
      double m = 1.23;
      gsl_interval x = {0.0, 2*M_PI};

      gsl_min_fminimizer * s 
	= gsl_min_fminimizer_alloc (gsl_min_fminimizer_brent);

      gsl_min_fminimizer_set (s, &line.f, m, x);

      do
	{
	iterations++;
	status = gsl_min_fminimizer_iterate (s);
	
	m = gsl_min_fminimizer_minimum (s);
	x = gsl_min_fminimizer_interval (s);
	
	status = gsl_min_test_interval (x, 0.001, 0.0);
	
	if (status == GSL_SUCCESS)
	  printf ("Converged:\n");
	
	printf ("%5d [%.7f, %.7f] %.7f %.7f %+.7f\n",
		iterations, x.lower, x.upper, 
		m, s->f_minimum, x.upper - x.lower);
	}
      while (status == GSL_CONTINUE && iterations < 100);
    }
  }

      
  return 0;
}