#include <stdlib.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_errno.h>

double pi;

double spring(const gsl_vector *x, void *pa) {
  double x0=gsl_vector_get(x,0);
  double x1=gsl_vector_get(x,1);
  double x2=gsl_vector_get(x,2);

  double theta=atan2(x1,x0);
  double r=sqrt(x0*x0+x1*x1);
  double z=x2;
  while (z>pi) z-=2.0*pi;
  while (z<-pi) z+=2.0*pi;
  double tmz=theta-z;
  double rm1=r-1.0;
  double ret=exp(tmz*tmz+rm1*rm1)+fabs(x2/10.0);
  if (!finite(ret)) exit(-1);
  return ret;
}

int main(int argv, char *argc[]) {

  pi=acos(-1.0);

  int nvar=3;

  {
    const gsl_multimin_fminimizer_type *T=gsl_multimin_fminimizer_nmsimplex;
    gsl_multimin_fminimizer *s=0;
    gsl_multimin_function minex_func;

    gsl_vector *gx=gsl_vector_alloc(nvar);
    gsl_vector_set(gx,0,1.0);
    gsl_vector_set(gx,1,0.0);
    gsl_vector_set(gx,2,7.0*pi);
	
    int iter=0;
    int status;
    double size;
     
    gsl_vector *ss;
    ss=gsl_vector_alloc(nvar);
    gsl_vector_set_all(ss,1.0);
     
    minex_func.n=nvar;
    minex_func.f=&spring;
    minex_func.params=0;
      
    s=gsl_multimin_fminimizer_alloc(T,nvar);
    gsl_multimin_fminimizer_set(s,&minex_func,gx,ss);
      
    do {
      iter++;
      status=gsl_multimin_fminimizer_iterate(s);
        
      if (status) break;

      size=gsl_multimin_fminimizer_size(s);
      status=gsl_multimin_test_size(size,1.0e-4);

      if (iter%20==0) {
	printf("%d %6.5e %6.5e %6.5e\n",iter,
	       gsl_vector_get(s->x,0),
	       gsl_vector_get(s->x,1),
	       gsl_vector_get(s->x,2));
      }
	
    } while (status == GSL_CONTINUE && iter < 600);

    gsl_vector_free(ss);
    gsl_vector_free(gx);
    gsl_multimin_fminimizer_free(s);
  }

  printf("\n");

  {
    const gsl_multimin_fminimizer_type *T=gsl_multimin_fminimizer_nmsimplex2;
    gsl_multimin_fminimizer *s=0;
    gsl_multimin_function minex_func;

    gsl_vector *gx=gsl_vector_alloc(nvar);
    gsl_vector_set(gx,0,1.0);
    gsl_vector_set(gx,1,0.0);
    gsl_vector_set(gx,2,7.0*pi);
	
    int iter=0;
    int status;
    double size;
     
    gsl_vector *ss;
    ss=gsl_vector_alloc(nvar);
    gsl_vector_set_all(ss,1.0);
     
    minex_func.n=nvar;
    minex_func.f=&spring;
    minex_func.params=0;
      
    s=gsl_multimin_fminimizer_alloc(T,nvar);
    gsl_multimin_fminimizer_set(s,&minex_func,gx,ss);
      
    do {
      iter++;
      status=gsl_multimin_fminimizer_iterate(s);
        
      if (status) break;

      size=gsl_multimin_fminimizer_size(s);
      status=gsl_multimin_test_size(size,1.0e-4);

      if (iter%20==0) {
	printf("%d %6.5e %6.5e %6.5e\n",iter,
	       gsl_vector_get(s->x,0),
	       gsl_vector_get(s->x,1),
	       gsl_vector_get(s->x,2));
      }
	
    } while (status == GSL_CONTINUE && iter < 600);

    gsl_vector_free(ss);
    gsl_vector_free(gx);
    gsl_multimin_fminimizer_free(s);
  }

  return 0;
}