#include <iostream>
#include <string.h>
#include <fstream>
#include <math.h>
#include <time.h>
#include <cstdlib>
#include <ctime>
#include <conio.h>
#include <stdlib.h>
#include <stdio.h>
//R_A_N_2_P_A_R_A_M_E_T_E_R_S___________________________________________________
#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0 - EPS)
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
#define Pi 3.14159265359
using namespace std;
#define T 1024*4
#define J 1000///discretization for theoretical frequency calculations
#define Q 1//Discretization for parameter scanning
#define P 1//Trials
#define dt 0.1
#define Dt 0.001
#define N 100
#define alpha 1.0
#define gain 300
#define h 0.0
#define sparseness 0.
#define g -10
#define delay 25
#define R_c -1.375
#define D_max 0.1
// PROBLEM alpha^-1 *dx/dt = -x+R_1*x(t-/tau)
double RealEigenproblem(double rel, double iml,double R_1, double tau);
double Heaviside(double u);
double F_corrected(double u, double A, double w);
double F_Axel(double u, double A, double w);
double ImaginaryEigenproblem(double rel, double iml,double R_1, double tau);
float ran2(long *idum);//random number generator - initial cond.
void shuffle(double Array[], int size);//shuffle an array with double precision entries
void four1(double data[], unsigned long nn, int isign);// FFT NR routine
double signum(double u);
double f(double u);
double f_prime(double u);
double f_2prime(double u);
double f_3prime(double u);
double f_4prime(double u);
double f_5prime(double u);
double f_7prime(double u);
double f_9prime(double u);
double f_11prime(double u);
double u[N][T];
double u_spike[N][T];
double u_spatial_ave[T];
double u_ave[T];
double u_lin[T];
double u_approx[T];
double u_ave_spike[T];
double u_mean[T];
double u_fixed_point[T];
double nu[T];
double Dummy_1[T];
double Dummy_2[T];
double Dummy_3[T];
double Dummy_4[T];
double v[N][T];
double F[N][T];
double F2[N][T];
double F3[N][T];
double X[N][T];
double XI[N][T];
double W[N][N];
double PSD_1[(int)T];//Power Spectral Density
double PSD_2[(int)T];//Power Spectral Density
double PSD_3[(int)T];//Power Spectral Density
double PSD_4[(int)T];//Power Spectral Density
double Freq[T/2];//Array of output frequencies for PSD display
double D[Q];
double I_o[Q];
double freq[Q];
double PERIODIC[N][T];
double Peak_frequency_1[Q];
double Peak_frequency_2[Q];
double Peak_frequency_3[Q];
double Peak_frequency_4[Q];
double Theoretical_frequency[Q];
double Susceptibility[Q];
double Susceptibility_theo[Q];
double Fixed_point[Q];
double Fixed_point_theo[Q];
int main()
{
double tau= delay*dt;
double mean_w=0;
for (int n=0;n<N;n++)
{
for (int l=0;l<N;l++)
{
long seed_w=(long) n*l+l+n+34;
long seed=(long) 32+n*l+n+l+12;
if(ran2(&seed)<sparseness)
{
W[n][l]=0;
// W[l][n]=W[n][l];
}
//else{W[n][l] = -g*ran2(&seed_w); /*W[l][n]=W[n][l];*/ }
else{W[n][l]=g+2*fabs(g)*sqrt(-2*log(ran2(&seed_w)))*cos(2*3.14159*ran2(&seed));}
mean_w=mean_w+1/(((double) N)*((double)N))*W[n][l];
}
}
// cout<<mean_w<<endl;
for (int q=0;q<Q;q++)
{
cout<<"WARNING: includes periodic forcing!"<<endl;
D[q] =0.0;//+ D_max*q/((double)Q-1);
I_o[q]=1.0;//+1*q/((double)Q);
freq[q]=100;
cout<<freq[q]<<endl;
double omega = freq[q];//(2*Pi)*freq[q]*Dt/dt;
cout<<I_o[q]/freq[q]/*I_o[q]*I_o[q]*Pi/(((omega*omega)+1)*omega)*/<<endl;
for (int p=0;p<P;p++)
{
double var_nu=0;
for (int t=0;t<T;t++)
{
for (int n=0;n<N;n++)
{
long d=rand()%105;
long noisyseed1=(long)21*t+d*n+q*p+90*p+45*q;
long noisyseed2=(long)69*t+11*n+d*t+q*p+901*p+415*q;
XI[n][t] = sqrt(-2*log(ran2(&noisyseed1)))*cos(2*3.14159*ran2(&noisyseed2));
PERIODIC[n][t] = I_o[q]*sin(2*Pi*freq[q]*t*Dt);
}
nu[t+1]=nu[t]+dt*((-nu[t])+I_o[q]*PERIODIC[0][t])+sqrt(2*D[q]*dt)*XI[0][t];
var_nu=var_nu+1/((double)T)*(nu[t+1]*nu[t+1]);
}
// cout<<"var nu:"<<var_nu<<endl;
for (int t=tau;t<T;t++)
{
for (int i=0;i<N;i++)
{
u[i][t+1]= u[i][t]+dt*(-u[i][t]+F[i][t-delay]+PERIODIC[i][t]);
v[i][t+1]= v[i][t]+dt*(-v[i][t]+F2[i][t-delay]);
u_spike[i][t+1]=u_spike[i][t]+dt*(-u_spike[i][t]+F3[i][t-delay]+PERIODIC[i][t]);
F[i][t]=0;
F2[i][t]=0;
F3[i][t]=0;
for (int j=0;j<N;j++)
{
F[i][t] = F[i][t] + 1/((double)N)*W[i][j]*f(u[j][t]);
}
for (int j=0;j<N;j++)
{
F2[i][t] = F2[i][t] + 1/((double)N)*W[i][j]*f(v[j][t]);
}
for (int j=0;j<N;j++)
{
F3[i][t] = F3[i][t] + 1/((double)N)*W[i][j]*X[j][t];
}
long seed =rand()%101;
long seed_1 = (long) 3+seed*i+11*q+45*p*i+i+t;
double p =ran2(&seed_1);
double p_fire = 1-exp(-f(u_spike[i][t+1])*dt);
if (p<p_fire)//spike occurs
{
X[i][t+1] =1/dt;
}
else//no spike
{
X[i][t+1]=0;
}
u_spatial_ave[t+1] = u_spatial_ave[t+1]+1/((double)N)*u[i][t+1];
}
u_ave_spike[t]=0;
for (int i=0;i<N;i++)
{
u_ave_spike[t]=u_ave_spike[t]+1/((double)N)*u_spike[i][t+1];
}
double A1 = 0;
double B = 0;
for (int s=0;s<T;s++)
{
A1 = A1 + 1/((double) T)*f(u_ave[t-delay]+nu[s]);
B = B + 1/((double) T)*f(u_fixed_point[t]+nu[s]);
}
//cout<<C<<endl;
u_ave[t+1]= u_ave[t]+dt*(-u_ave[t]+mean_w*A1);
u_mean[t+1]= u_mean[t]+dt*(-u_mean[t]+mean_w*(F_Axel(u_mean[t-delay], I_o[q],freq[q]/(2*Pi)))+1*PERIODIC[0][t]);
//u_mean[t+1]= u_mean[t]+dt*(-u_mean[t]+mean_w*(f(u_mean[t-delay]))+PERIODIC[0][t]);
u_fixed_point[t+1] = u_fixed_point[t]+dt*(-u_fixed_point[t]+mean_w*B);
}
double R;
double c_min=-I_o[q]/freq[q]+0.00001;
int Nc =1000;
double c_max=I_o[q]/freq[q]-0.00001;
double dc = fabs(c_min-c_max)/((double) Nc);
double C=0;
double E=0;
double F=0;
double G=0;
for (int k=0;k<Nc;k++)
{
double c = c_min+dc*k;
C = C + dc*mean_w*f_prime(u_fixed_point[T-10]+c)*1/sqrt((I_o[q]/freq[q]*I_o[q]/freq[q])-c*c)*1/Pi;
}
R=C;
Fixed_point[q]=Fixed_point[q]+1/((double)P)*u_fixed_point[T-10];
Fixed_point_theo[q]=0;//-sqrt(Pi)*sqrt(D[q])*mean_w/(sqrt(2)*g-2*sqrt(Pi)*sqrt(D[q]));
// Susceptibility[q]=Susceptibility[q]+1/((double)P)*R;
Susceptibility[q]=Susceptibility[q]+1/((double)P)*(C);
// Susceptibility_theo[q]=mean_w*(Heaviside(Fixed_point[q]+l)/(sqrt(0.1e-1-(Fixed_point[q]*Fixed_point[q])+l^2)*Pi)-Heaviside(mu-l)/(sqrt(0.1e-1-mu^2+l^2)*Pi))
// cout<<"R:"<<Susceptibility[q]<<endl;
// cout<<"fp:"<<Fixed_point[q]<<endl;
for (int t=tau;t<T;t++)
{
u_lin[t+1]=u_lin[t]+dt*(-u_lin[t]+R*u_lin[t-delay]);
}
//Compute mean EEG
//Compute mean EEG
double mean_1=0;
double mean_2=0;
double mean_3=0;
double mean_4=0;
for (int t=0;t<(int)T;t++)
{
mean_1 = mean_1+1/((double) T)*u[(int)N/2][t];
mean_2 = mean_2+1/((double) T)*v[(int)N/2][t];
mean_3 = mean_3+1/((double) T)*u_mean[t];
mean_4 = mean_4+1/((double) T)*u_ave_spike[t];
}
//Compute mean corrected EEG signal
for (int t=0;t<(int)T;t++)
{
Dummy_1[t] = u[(int)N/2][t]-mean_1;
Dummy_2[t] = v[(int)N/2][t]-mean_2;
Dummy_3[t] = u_mean[t]-mean_3;
Dummy_4[t] = u_ave_spike[t]-mean_4;
}
//fourier transforms and power spectral denstities
for (int k=0;k<(int)T;k++)
{
PSD_1[k] = 0;
PSD_2[k] = 0;
PSD_3[k] = 0;
PSD_4[k] = 0;
}
for (int k=0;k<(int)T;k++)
{
PSD_1[k] = Dummy_1[k];
PSD_2[k] = Dummy_2[k];
PSD_3[k] = Dummy_3[k];
PSD_4[k] = Dummy_4[k];
}
unsigned long nn=T/2;
four1(PSD_1-1, nn,1);
for (int k=0;k<T/2;k++)
{
PSD_1[k] = 1/((double)T*T)*(fabs(PSD_1[k])*fabs(PSD_1[k])+fabs(PSD_1[(int)T-k])*fabs(PSD_1[(int)T-k]));
}
four1(PSD_2-1, nn,1);
for (int k=0;k<T/2;k++)
{
PSD_2[k] = 1/((double)T*T)*(fabs(PSD_2[k])*fabs(PSD_2[k])+fabs(PSD_2[(int)T-k])*fabs(PSD_2[(int)T-k]));
}
four1(PSD_3-1, nn,1);
for (int k=0;k<T/2;k++)
{
PSD_3[k] = 1/((double)T*T)*(fabs(PSD_3[k])*fabs(PSD_3[k])+fabs(PSD_3[(int)T-k])*fabs(PSD_3[(int)T-k]));
}
four1(PSD_4-1, nn,1);
for (int k=0;k<T/2;k++)
{
PSD_4[k] = 1/((double)T*T)*(fabs(PSD_4[k])*fabs(PSD_4[k])+fabs(PSD_4[(int)T-k])*fabs(PSD_4[(int)T-k]));
}
double max_f1=0;
double max_power1=0;
double max_f2=0;
double max_power2=0;
double max_f3=0;
double max_power3=0;
double max_f4=0;
double max_power4=0;
for(int k =2*2*(T)*Dt; k<2*50*(T)*Dt;k++)
{
if ( PSD_1[k] >max_power1)
{
max_power1 = PSD_1[k];
max_f1=k*1/((double) 2* (T)*Dt);
}
}
for(int k =2*2*(T)*Dt; k<2*50*(T)*Dt;k++)
{
if ( PSD_2[k] >max_power2)
{
max_power2 = PSD_2[k];
max_f2=k*1/((double) 2* (T)*Dt);
}
}
for(int k =2*2*(T)*Dt; k<2*50*(T)*Dt;k++)
{
if ( PSD_3[k] >max_power3)
{
max_power3 = PSD_3[k];
max_f3=k*1/((double) 2* (T)*Dt);
}
}
for(int k =2*2*(T)*Dt; k<2*50*(T)*Dt;k++)
{
if ( PSD_4[k] >max_power4)
{
max_power4 = PSD_4[k];
max_f4=k*1/((double) 2* (T)*Dt);
}
}
Peak_frequency_1[q] = Peak_frequency_1[q]+1/((double) P)*max_f1;
Peak_frequency_2[q] = Peak_frequency_2[q]+1/((double) P)*max_f2;
Peak_frequency_3[q] = Peak_frequency_3[q]+1/((double) P)*max_f3;
Peak_frequency_4[q] = Peak_frequency_4[q]+1/((double) P)*max_f4;
Theoretical_frequency[q] = -(dt/Dt)/(2*Pi)*(-Pi+acos((-1)/R_c))/(delay*dt);;//-2*Pi*R*sqrt(1-(-1/R)*(-1/R));//
cout<<Peak_frequency_1[q]<<" "<<Peak_frequency_2[q]<<endl;
}//p loop
}//q loop
double H=0;
int No_of_points=100;
double U_max=0.5;
double U_min=-0.5;
double U[No_of_points];
double F[No_of_points];
double F_approx[No_of_points];
int Nc =100000;
double Input_Level=0.1;
double fre =2*Pi*50/100;
double c_min=-Input_Level/fre+0.001*Input_Level;
double c_max=Input_Level/fre-0.001*Input_Level;
double dc = fabs(c_min-c_max)/((double) Nc);
for (int k1=0;k1<100;k1++)
{
U[k1] = U_min+fabs(U_max-U_min)/((double)No_of_points)*k1;
H=0;
for (int k2=0;k2<Nc;k2++)
{
double c = c_min+dc*k2;
// cout<<c<<endl;
H = H + dc*f(U[k1]+c)*1/(Pi*sqrt(Input_Level/fre*Input_Level/fre-c*c));
}
cout<<H<<endl;
F[k1] = H;
//F_approx[k1]=0.5+(0.5)*erf((0.5)*sqrt(2)*U[k1]/sqrt(Noise_Level));
}
for(int k=0;k<(int) T/2;k++)
{
Freq[k] = k*1/((double) 2* (T)*Dt);
//cout<<k*1/((double) 2* (T/2)*Dt)<<endl;
}
ofstream outfile;
outfile.open("SS - PSD.txt", ios::out);
for(int k=0;k<(int) 30*(2*T)*Dt;k++)
{
outfile<<Freq[k]<<" "<<PSD_1[k]<<" "<<PSD_2[k]<<endl;
}
outfile.close();
outfile.open("SS - Connectivity.txt", ios::out);
for(int i=0;i<N;i++)
{
for (int j=0;j<N;j++)
{
outfile<<i<<" "<<j<<" "<<W[i][j]<<endl;
}
}
outfile.close();
outfile.open("SS - Noise induced Linearization.txt", ios::out);
for(int k=0;k<No_of_points;k++)
{
outfile<<U[k]<<" "<<F[k]<<" "<<F_approx[k]<<endl;
}
outfile.close();
outfile.open("SS - Sliced view.txt", ios::out);
for(int i=0;i<T;i++)
{
outfile<<i*Dt<<" "<<u[(int)N/2][i]<<" "<<v[(int)N/2][i]<<" "<< u_spatial_ave[i]<<" "<<u_mean[i]<<" "<<u_ave_spike[i]<<" "<<I_o[0]*PERIODIC[(int)N/2][i]<<endl;
}
outfile.close();
outfile.open("SS - misceleanous view.txt", ios::out);
for(int i=0;i<T;i++)
{
outfile<<i*Dt<<" "<<u[(int)N/2][i]<<" "<<u[(int)N/2+15][i]<<" "<<u[(int)N/2+9][i]<<" "<<u[(int)N/2-15][i]<<" "<<u[(int)N/2+10][i]<<endl;
}
outfile.close();
outfile.open("SS - Spikes.txt", ios::out);
for(int i=0;i<T;i++)
{
for (int n=0;n<N;n++)
{
if (X[n][i])
{
outfile<<i*Dt<<" "<<n<<endl;
}
}
}
outfile.close();
outfile.open("SS - Peak frequency.txt", ios::out);
for(int q=0;q<Q;q++)
{
outfile<<I_o[q]/freq[q]<<" "<<Peak_frequency_1[q]<<" "<<Peak_frequency_2[q]<<" "<<Peak_frequency_3[q]<<" "<<Peak_frequency_4[q]<<" "<<Theoretical_frequency[q]<<endl;
}
outfile.close();
outfile.open("SS - Corrected F with Periodic Forcing.txt", ios::out);
for(int q=0;q<1000;q++)
{
double input=-0.5+1*q/((double)1000);
outfile<<input<<" "<<F_corrected(input, 0.01,2*Pi*50/(100))<<endl;
}
outfile.close();
outfile.open("SS - R.txt", ios::out);
for(int q=0;q<Q;q++)
{
outfile<<I_o[q]<<" "<<Susceptibility[q]<<" "<<Susceptibility_theo[q]<<endl;
}
outfile.close();
outfile.open("SS - Fixed Point.txt", ios::out);
for(int q=0;q<Q;q++)
{
outfile<<D[q]<<" "<<Fixed_point[q]<<" "<<Fixed_point_theo[q]<<endl;
}
outfile.close();
system("pause");
cout<<endl;
return 0;
}
double RealEigenproblem(double rel, double iml, double R_1, double tau)
{
double output=rel*(1/alpha)+1+R_1*exp(-rel*tau)*cos(iml*tau);
return output;
}
double ImaginaryEigenproblem(double rel, double iml,double R_1, double tau)
{
double output=iml*(1/alpha)-R_1*exp(-rel*tau)*sin(iml*tau);
return output;
}
double signum(double u)
{
double output;
if (u>0)
{
output=1;
}
else{output=-1;}
return output;
}
double f(double u)
{
double output;
output = 1/(1+exp(-gain*(u-h)));
//output = tanh(u);//1/(1+exp(-gain*(u-h)));
return output;
}
double f_prime(double u)
{
double output;
output = pow(0.1e1 + exp(-gain * (u - h)), -0.2e1) * gain * exp(-gain * (u - h));
//output = 1-tanh(u)*tanh(u);//pow(0.1e1 + exp(-gain * (u - h)), -0.2e1) * gain * exp(-gain * (u - h));
return output;
}
double f_2prime(double u)
{
double output;
output = 0.2e1 * pow(0.1e1 + exp(-gain * (u - h)), -0.3e1) * gain * gain * pow(exp(-gain * (u - h)), 0.2e1) - pow(0.1e1 + exp(-gain * (u - h)), -0.2e1) * gain * gain * exp(-gain * (u - h));
//output = -0.2e1 * tanh(u) * (0.1e1 - pow(tanh(u), 0.2e1)); //0.2e1 * pow(0.1e1 + exp(-gain * (u - h)), -0.3e1) * gain * gain * pow(exp(-gain * (u - h)), 0.2e1) - pow(0.1e1 + exp(-gain * (u - h)), -0.2e1) * gain * gain * exp(-gain * (u - h));
return output;
}
double f_3prime(double u)
{
double output;
output = 0.6e1 * pow(0.1e1 + exp(-gain * (u - h)), -0.4e1) * pow(gain, 0.3e1) * pow(exp(-gain * (u - h)), 0.3e1) - 0.6e1 * pow(0.1e1 + exp(-gain * (u - h)), -0.3e1) * pow(gain, 0.3e1) * pow(exp(-gain * (u - h)), 0.2e1) + pow(0.1e1 + exp(-gain * (u - h)), -0.2e1) * pow(gain, 0.3e1) * exp(-gain * (u - h));
//output = -0.2e1 * pow(0.1e1 - pow(tanh(u), 0.2e1), 0.2e1) + 0.4e1 * pow(tanh(u), 0.2e1) * (0.1e1 - pow(tanh(u), 0.2e1));// 0.6e1 * pow(0.1e1 + exp(-gain * (u - h)), -0.4e1) * pow(gain, 0.3e1) * pow(exp(-gain * (u - h)), 0.3e1) - 0.6e1 * pow(0.1e1 + exp(-gain * (u - h)), -0.3e1) * pow(gain, 0.3e1) * pow(exp(-gain * (u - h)), 0.2e1) + pow(0.1e1 + exp(-gain * (u - h)), -0.2e1) * pow(gain, 0.3e1) * exp(-gain * (u - h));
return output;
}
double f_5prime(double u)
{
double output;
output = 0.120e3 * pow(0.1e1 + exp(-gain * (u - h)), -0.6e1) * pow(gain, 0.5e1) * pow(exp(-gain * (u - h)), 0.5e1) - 0.240e3 * pow(0.1e1 + exp(-gain * (u - h)), -0.5e1) * pow(gain, 0.5e1) * pow(exp(-gain * (u - h)), 0.4e1) + 0.150e3 * pow(0.1e1 + exp(-gain * (u - h)), -0.4e1) * pow(gain, 0.5e1) * pow(exp(-gain * (u - h)), 0.3e1) - 0.30e2 * pow(0.1e1 + exp(-gain * (u - h)), -0.3e1) * pow(gain, 0.5e1) * pow(exp(-gain * (u - h)), 0.2e1) + pow(0.1e1 + exp(-gain * (u - h)), -0.2e1) * pow(gain, 0.5e1) * exp(-gain * (u - h));
return output;
}
double f_4prime(double u)
{
double output;
output = 0.24e2 * pow(0.1e1 + exp(-gain * (u - h)), -0.5e1) * pow(gain, 0.4e1) * pow(exp(-gain * (u - h)), 0.4e1) - 0.36e2 * pow(0.1e1 + exp(-gain * (u - h)), -0.4e1) * pow(gain, 0.4e1) * pow(exp(-gain * (u - h)), 0.3e1) + 0.14e2 * pow(0.1e1 + exp(-gain * (u - h)), -0.3e1) * pow(gain, 0.4e1) * pow(exp(-gain * (u - h)), 0.2e1) - pow(0.1e1 + exp(-gain * (u - h)), -0.2e1) * pow(gain, 0.4e1) * exp(-gain * (u - h));
return output;
}
double f_7prime(double u)
{
double output;
output = 0.5040e4 * pow(0.1e1 + exp(-gain * (u - h)), -0.8e1) * pow(gain, 0.7e1) * pow(exp(-gain * (u - h)), 0.7e1) - 0.15120e5 * pow(0.1e1 + exp(-gain * (u - h)), -0.7e1) * pow(gain, 0.7e1) * pow(exp(-gain * (u - h)), 0.6e1) + 0.16800e5 * pow(0.1e1 + exp(-gain * (u - h)), -0.6e1) * pow(gain, 0.7e1) * pow(exp(-gain * (u - h)), 0.5e1) - 0.8400e4 * pow(0.1e1 + exp(-gain * (u - h)), -0.5e1) * pow(gain, 0.7e1) * pow(exp(-gain * (u - h)), 0.4e1) + 0.1806e4 * pow(0.1e1 + exp(-gain * (u - h)), -0.4e1) * pow(gain, 0.7e1) * pow(exp(-gain * (u - h)), 0.3e1) - 0.126e3 * pow(0.1e1 + exp(-gain * (u - h)), -0.3e1) * pow(gain, 0.7e1) * pow(exp(-gain * (u - h)), 0.2e1) + pow(0.1e1 + exp(-gain * (u - h)), -0.2e1) * pow(gain, 0.7e1) * exp(-gain * (u - h));
return output;
}
double f_9prime(double u)
{
double output;
output = 0.362880e6 * pow(0.1e1 + exp(-gain * (u - h)), -0.10e2) * pow(gain, 0.9e1) * pow(exp(-gain * (u - h)), 0.9e1) - 0.1451520e7 * pow(0.1e1 + exp(-gain * (u - h)), -0.9e1) * pow(gain, 0.9e1) * pow(exp(-gain * (u - h)), 0.8e1) + 0.2328480e7 * pow(0.1e1 + exp(-gain * (u - h)), -0.8e1) * pow(gain, 0.9e1) * pow(exp(-gain * (u - h)), 0.7e1) - 0.1905120e7 * pow(0.1e1 + exp(-gain * (u - h)), -0.7e1) * pow(gain, 0.9e1) * pow(exp(-gain * (u - h)), 0.6e1) + 0.834120e6 * pow(0.1e1 + exp(-gain * (u - h)), -0.6e1) * pow(gain, 0.9e1) * pow(exp(-gain * (u - h)), 0.5e1) - 0.186480e6 * pow(0.1e1 + exp(-gain * (u - h)), -0.5e1) * pow(gain, 0.9e1) * pow(exp(-gain * (u - h)), 0.4e1) + 0.18150e5 * pow(0.1e1 + exp(-gain * (u - h)), -0.4e1) * pow(gain, 0.9e1) * pow(exp(-gain * (u - h)), 0.3e1) - 0.510e3 * pow(0.1e1 + exp(-gain * (u - h)), -0.3e1) * pow(gain, 0.9e1) * pow(exp(-gain * (u - h)), 0.2e1) + pow(0.1e1 + exp(-gain * (u - h)), -0.2e1) * pow(gain, 0.9e1) * exp(-gain * (u - h));
return output;
}
double f_11prime(double u)
{
double output;
output = 0.39916800e8 * pow(0.1e1 + exp(-gain * (u - h)), -0.12e2) * pow(gain, 0.11e2) * pow(exp(-gain * (u - h)), 0.11e2) - 0.199584000e9 * pow(0.1e1 + exp(-gain * (u - h)), -0.11e2) * pow(gain, 0.11e2) * pow(exp(-gain * (u - h)), 0.10e2) + 0.419126400e9 * pow(0.1e1 + exp(-gain * (u - h)), -0.10e2) * pow(gain, 0.11e2) * pow(exp(-gain * (u - h)), 0.9e1) - 0.479001600e9 * pow(0.1e1 + exp(-gain * (u - h)), -0.9e1) * pow(gain, 0.11e2) * pow(exp(-gain * (u - h)), 0.8e1) + 0.322494480e9 * pow(0.1e1 + exp(-gain * (u - h)), -0.8e1) * pow(gain, 0.11e2) * pow(exp(-gain * (u - h)), 0.7e1) - 0.129230640e9 * pow(0.1e1 + exp(-gain * (u - h)), -0.7e1) * pow(gain, 0.11e2) * pow(exp(-gain * (u - h)), 0.6e1) + 0.29607600e8 * pow(0.1e1 + exp(-gain * (u - h)), -0.6e1) * pow(gain, 0.11e2) * pow(exp(-gain * (u - h)), 0.5e1) - 0.3498000e7 * pow(0.1e1 + exp(-gain * (u - h)), -0.5e1) * pow(gain, 0.11e2) * pow(exp(-gain * (u - h)), 0.4e1) + 0.171006e6 * pow(0.1e1 + exp(-gain * (u - h)), -0.4e1) * pow(gain, 0.11e2) * pow(exp(-gain * (u - h)), 0.3e1) - 0.2046e4 * pow(0.1e1 + exp(-gain * (u - h)), -0.3e1) * pow(gain, 0.11e2) * pow(exp(-gain * (u - h)), 0.2e1) + pow(0.1e1 + exp(-gain * (u - h)), -0.2e1) * pow(gain, 0.11e2) * exp(-gain * (u - h));
return output;
}
float ran2(long *idum)
{
int j;
long k;
static long idum2 = 123456789;
static long iy = 0;
static long iv[NTAB];
float temp;
if (*idum <= 0) { /* initialize */
if (-(*idum) < 1) /* prevent idum == 0 */
*idum = 1;
else
*idum = -(*idum); /* make idum positive */
idum2 = (*idum);
for (j = NTAB + 7; j >= 0; j--) { /* load the shuffle table */
k = (*idum) / IQ1;
*idum = IA1 * (*idum - k*IQ1) - k*IR1;
if (*idum < 0)
*idum += IM1;
if (j < NTAB)
iv[j] = *idum;
}
iy = iv[0];
}
k = (*idum) / IQ1;
*idum = IA1 * (*idum - k*IQ1) - k*IR1;
if (*idum < 0)
*idum += IM1;
k = idum2/IQ2;
idum2 = IA2 * (idum2 - k*IQ2) - k*IR2;
if (idum2 < 0)
idum2 += IM2;
j = iy / NDIV;
iy = iv[j] - idum2;
iv[j] = *idum;
if (iy < 1)
iy += IMM1;
if ((temp = AM * iy) > RNMX)
return RNMX; /* avoid endpoint */
else
return temp;
}
/******************************************************************************/
void four1(double data[], unsigned long nn, int isign)
/*******************************************************************************
Replaces data[1..2*nn] by its discrete Fourier transform, if isign is input as
1; or replaces data[1..2*nn] by nn times its inverse discrete Fourier transform,
if isign is input as -1. data is a complex array of length nn or, equivalently,
a real array of length 2*nn. nn MUST be an integer power of 2 (this is not
checked for!).
*******************************************************************************/
{
unsigned long n,mmax,m,j,istep,i;
double wtemp,wr,wpr,wpi,wi,theta;
double tempr,tempi;
n=nn << 1;
j=1;
for (i=1;i<n;i+=2) { /* This is the bit-reversal section of the routine. */
if (j > i) {
SWAP(data[j],data[i]); /* Exchange the two complex numbers. */
SWAP(data[j+1],data[i+1]);
}
m=nn;
while (m >= 2 && j > m) {
j -= m;
m >>= 1;
}
j += m;
}
mmax=2;
while (n > mmax) { /* Outer loop executed log2 nn times. */
istep=mmax << 1;
theta=isign*(6.28318530717959/mmax); /* Initialize the trigonometric recurrence. */
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0;
wi=0.0;
for (m=1;m<mmax;m+=2) { /* Here are the two nested inner loops. */
for (i=m;i<=n;i+=istep) {
j=i+mmax; /* This is the Danielson-Lanczos formula. */
tempr=wr*data[j]-wi*data[j+1];
tempi=wr*data[j+1]+wi*data[j];
data[j]=data[i]-tempr;
data[j+1]=data[i+1]-tempi;
data[i] += tempr;
data[i+1] += tempi;
}
wr=(wtemp=wr)*wpr-wi*wpi+wr; /* Trigonometric recurrence. */
wi=wi*wpr+wtemp*wpi+wi;
}
mmax=istep;
}
}
double Heaviside(double u)
{
double output;
if (u<0)
{
output=0;
}
else
{
output=1;
}
return output;
}
double F_corrected(double u, double A, double w)
{
double output;
if (u>-A/w&&u<A/w)
{
output=(0.5)+asin(w*u/A)/Pi;
}
else
{
if (u<0)
{
output=0;
}
else
{
output=1;
}
}
return output;
}
double F_Axel(double u, double A, double w)
{
double output;
double mu=A/w;
if (u>-A/w&&u<A/w)
{
output=0.5*((2/Pi)*asin(u/mu)+1);
}
else
{
if (u<-mu)
{
output=0;
}
else
{
output=1;
}
}
return output;
}