/////gradient.c
/////傾きが最大となる点を検出するプログラム
/////2013/12/25 by Takashi Kojima

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

#define N 1000

int main(){

  int    i;
  double x[N];
  double y[N];
  double grad[N];
  double ave_grad[N];
  double max_grad, max_x;
  double min_grad, min_x;
  double abs_max_grad, abs_max_x;
  int    data_length=0;

  //観測データの読み込み
  FILE *fp;
  char *filename1 = "s0.1_d0.5.dat";

  fp = fopen(filename1, "r");
  if( fp == NULL ){
    printf("Can't open file %s\n", filename1);
    exit(1);
  }
  i=0;
  while(fscanf(fp, "%lf %lf", &x[i], &y[i]) != EOF){
    data_length++;
    i++;
  }

  printf("Finished reaing.\n");
  fclose(fp);
  

  //勾配計算とその周囲4つによる平均化
  for(i=0 ; i<data_length ; i++){
    grad[i]=(y[i+1]-y[i])/(x[i+1]-x[i]);
  }

  for(i=2 ; i<data_length-1 ; i++){
    ave_grad[i] = (grad[i+1] + grad[i] + grad[i-1] + grad[i-2])/4.0;
  }


  //データ書き出し
  //FILE *fp;
  char *filename2 = "grad_s0.1_d0.5.dat";

  fp = fopen(filename2, "w");
  if( fp==NULL ){
    printf("Can't open file %s\n", filename2);
    exit(1);
  }

  for(i=2 ; i<data_length-1 ; i++){
    //fprintf(fp, "%lf   %lf  %lf\n" ,x[i], ave_grad[i], grad[i]);
    fprintf(fp, "%lf   %lf \n" ,x[i], ave_grad[i]);
  }

  printf("Finished saving.\n");
  fclose(fp);


  //微係数の①最大値、②最小値、③絶対最大値の検出
  max_grad=0.0;
  min_grad=0.0;

  for (i=2; i<data_length-1 ; i++){
    if(grad[i]>=max_grad){
      max_grad = grad[i];
      max_x = x[i];
    }
    if(grad[i]<=min_grad){
      min_grad = grad[i];
      min_x = x[i];
    }
  }

  if( abs(max_grad) < abs(min_grad) ){
    abs_max_grad = min_grad;
    abs_max_x = max_x; 
  }else{
    abs_max_grad = max_grad;
    abs_max_x = min_x; 
  }
  
  printf("grad_max=%lf (x=%lf)\n", max_grad, max_x);
  printf("grad_min=%lf (x=%lf)\n", min_grad, min_x);
  printf("abs_grad_max=%lf (x=%lf)\n", abs_max_grad, abs_max_x);  

  return 0;
}