/*
 * bisection2.c -- 二分法 (bisection method) で方程式 f(x)=0 を解く
 *   コンパイル: gcc -o bisection2 bisection2.c -lm
 *   実行: ./bisection2
 *   入力は r=1 j=1 など
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double r;
double f(double x)
{
  return tan(x)-2*r*x/(x*x-r*r);
}

int main()
{
  int i,j, maxitr = 100;
  double a, b, c, eps;
  double fa, fb, fc,pai2;

  printf("r=");
  scanf("%lf",&r);

  printf("j=");
  scanf("%d",&j);

  if(j<1){
    printf("Please Input j>=1 int.\n");
    exit(0);
  }

  //pai2=pai/2
  pai2=2*atan(1.0);

  //auto
  if(r<pai2){
    if(j==1){
      a=r+0.0001;
      b=pai2-0.0001;
    }
    if(j>1){
      a=2*pai2*(j-1);
      b=2*pai2*(j-1)+pai2-0.000001;
    }
  }
  if(r>pai2){
    if(j==1){
      a=pai2+0.0001;
      b=pai2*2-0.0001;
    }
    if(j>1){
      a=j*2*pai2-pai2+0.0001;
      b=j*pai2*2-0.0001;
    }
  }
  if(r==pai2){
    if(j==1){
      a=pai2-0.0001;
      b=pai2+0.0001;
    }
    if(j>1){
      a=2*pai2*(j-1);
      b=2*pai2*(j-1)+pai2-0.00001;
    }
  }
  if(r<=0){
    printf("Please Input r>0.\n");
    exit(0);
  }

  eps=1e-15;

  fa = f(a);
  fb = f(b);
  if (fa * fb > 0.0) {
    printf(" f(α) f(β) > 0 なのであきらめます。\n");
    exit(0);
  }
  else {
    printf("f(%20.15f)=%9.2e, f(%20.15f)=%9.2e\n", a, fa, b, fb);
    for (i = 0; i < maxitr; i++) {
      c = (a + b) / 2;
      fc = f(c);
      if (fc == 0.0)
	break;
      else if (fa * fc <= 0.0) {
	/* 左側 [a,c] に根がある */
	b = c;
	fb = fc;
      }
      else {
	/* 左側 [a,c] には根がないかもしれない。[c,b] にあるはず */
	a = c;
	fa = fc;
      }
      printf ("f(%20.15f)=%9.2e, f(%20.15f)=%9.2e\n", a, fa, b, fb);
      if ((b - a) <= eps)
	break;
    }
    printf ("二分法による近似解 x = %20.15f\n", c);
    printf("lamda_%d=%20.15f\n",j,c*c);
  }

  //gnuplot

  FILE *hoge;
  hoge=popen("gnuplot -persist","w");
  fprintf(hoge,"set xrange[0:1]\n");
  fprintf(hoge,"plot cos(%f*x)+%f*sin(%f*x)/%f \n",c,r,c,c);

  pclose(hoge);

  return 0;
}