/*
  k_integration

   lin_k_integ.c

   2003.11.14 H.Akitaya

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

#define DELTAS 1000.0

float calc_dkds ( float knc, float k, float s );
float calc_knc_fromfunc ( float s );
double calc_dkds_strict ( double s, double k );

main(int argc, char *argv[] )
{
  FILE *f_knc;
  float s_from[200], s_to[200], s_mid[200], knc[200], sig_knc[200],
    err_knc[200];
  float npix, rej_npix;
  float sinteg_from, sinteg_to, k_init;
  float dkds, delta_s, knc_now;
  double s, k;
  double rk1, rk2, rk3, rk4, h;
  int n_data;
  int flag_found, sign;

  long j;
  long i, total_step;
  
  if ( argc != 6 ){
    printf( "useage : lin_kinteg (knc_filename) (from) (to) (k_nc_init) (delta_s)\n" );
    exit(1);
  }
  if( ( f_knc = fopen(argv[1],"r") ) == NULL ){
    printf( "input k_nc file is not found !!!\n" );
    exit(1);
  }
  sinteg_from = atof( argv[2] );
  sinteg_to = atof( argv[3] );
  k_init = atof( argv[4] );
  delta_s = atof( argv[5] );
  /*delta_s = DELTAS;*/

  /*  printf("%f %f %f\n", sinteg_from, sinteg_to, k_init );*/

  n_data = 0 ;
  while( fscanf( f_knc, "%e%e%e%e%e%e%e%e",
		 &s_from[n_data], &s_to[n_data], &s_mid[n_data],
		 &npix, &knc[n_data], &sig_knc[n_data],
		 &err_knc[n_data], &rej_npix )!= EOF){
    /*printf( "%f %f %f\n", s_to[n_data], s_from[n_data], s_mid[n_data]);*/
    n_data++;
  }

  if( n_data == 0){
    printf("file has no data !!\n" );
    exit(1);
  }

  /* integration */

  k = k_init;
  s = sinteg_from;

  sign=1;
  total_step = abs( sinteg_to - sinteg_from ) / delta_s; 
  if ( (sinteg_to - sinteg_from ) < 0 ) {
    delta_s *= -1;
    sign*=-1;
  }


   h = delta_s;

   printf("\n");
   
   j = 0;
   /*  for(i=0; i < total_step; i++ ){*/
   for( s = sinteg_from; ( ( s - sinteg_to )*sign <= 0 ); s+=delta_s ){
			  
    /* get knc */
    
    /*    flag_found = 0;
    for( j=0; j < n_data; j++){
      if( s >= s_from[j] && s < s_to[j] ){
	knc_now = knc[j];
	flag_found = 1;
      }
    }
    if( flag_found != 1 ){
      printf("data not found !!\n");
      exit(1);
      }*/
    /* Runge-Kutta Method (4th order) */
    rk1 = h * calc_dkds_strict( s, k );
    rk2 = h * calc_dkds_strict( s + h / 2.0 , k + rk1 / 2.0 );
    rk3 = h * calc_dkds_strict( s + h / 2.0 , k + rk2 / 2.0 );
    rk4 = h * calc_dkds_strict( s + h , k + rk3 );
        
    k = k + ( rk1 + 2.0*rk2 + 2.0*rk3 + rk4 ) / 6.0;
    s = s + h;

    j++;
    if( j%2000 == 0 ){
      printf( "%f %f %e %f\n", s, knc_now, dkds, k );
    j=0;
    }
  }
  
}


float calc_knc_fromfunc ( float s ){
  float knc;
  float a, sa, c, sc, d, sd, e, se, km, f, sf;
  km = 0.455;
  a = -9.813e-7;
  sa = 31738.9;
  c = -3.41831e-6;
  sc = 4408.96232;
  d = -0.000189062;
  sd = 420.517079;
  e = 0.000300422;
  se = 204.79792;
  f = -0.000205688;
  sf = 7.53919;
  
  if ( s >= sa ){
    knc = a*( s - sa ) + km;
  }else if ( s>= sc ){
    knc = km;
  }else if ( s >= 347.07669 ){
    knc = c*( s - sc ) + km;
  }else if ( s >= 288.11892 ){
    knc = d*( s - sd ) + km;
  }else if ( s > 124.63006 ){
    knc = e*( s - se ) + km;
  }else {
    knc = f*( s - sf ) + km;
  }
  return( knc );
}

float calc_dkds ( float knc, float k, float s )
{
  float dkds;
  dkds = ( knc - k ) / s;
  return( dkds );
}

double calc_dkds_strict ( double s, double k )
{
  double dkds;
  dkds = ( sqrt( k * calc_knc_fromfunc ( s ) ) - k ) / s;
  return( dkds );
}