// Programm harmonischerOszillator simuliert einen harmonischen Oszillator 
// x'' = - w^2 x, w = 2 PI/T, der in jedem Sampling-Schritt zufaellig
// gestoert wird.
//
// 13.10.03, Thomas Maiwald

#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <math.h>
#include "nr.h"
#include "nrutil.h"

// Unterschiedlicher Name des Time-Headers in Windows und Unix:
#ifdef WIN32
	#include <sys/timeb.h>
#else
	#include <sys/time.h>
#endif

#define PI 3.141592653589793

int main(int argc, char** argv)
{
	// Lese Uebergabeparameter:
	if(argc != 3)
    {
        printf("Aufruf: harmonischerOszillator T N\n");
        exit(1); 
    }
	double      T = atof(argv[1]);
    const int   N = atoi(argv[2]);
    
    // Programm beenden, falls N zu klein ist 
    if(N < 2)
    {
        printf("N muss groesser als 2 sein.\n");
        exit(1); 
    }

    // Aktuelle Millisekunden dienen als 'seed' fuer den Zufallszahlengenerator:
	long millisekunden = 0;
	#ifdef WIN32
		struct __timeb64 timebuffer;
		_ftime64( &timebuffer );
        millisekunden = timebuffer.millitm;
	#else 
		struct timeval tp;
		gettimeofday(&tp, NULL);
		millisekunden = tp.tv_usec;
	#endif
    long idum = - millisekunden;

	// Simuliere harmonischen Oszillator
	double a1   = 2 * cos(2*PI/T);
	double a2   = -1;
    double* x   = new double[N];
    x[0]        = 1.8;  // Startwerte
    x[1]        = 2; 
    cout << x[0] << endl;
    cout << x[1] << endl; 
    for(int i=2; i<N; i++)
    {
        x[i] = a1 * x[i-1] + a2 * x[i-2];
        //cout<<x[i]<<endl;
    }

    // Nun mit zufaelligen Stoerungen aus N(0,0.01)
    double* xMitStoerung = new double[N];
    xMitStoerung[0]      = 1.8;  // Startwerte
    xMitStoerung[1]      = 2; 
    cout << xMitStoerung[0] << endl;
    cout << xMitStoerung[1] << endl; 
    for(int i=2; i<N; i++)
    {
        xMitStoerung[i] = a1  * xMitStoerung[i-1] 
                        + a2  * xMitStoerung[i-2] 
                        + 0.1 * gasdev(&idum);
        //cout<<xMitStoerung[i]<<endl;
    }

    // Schaetzen der Varianz mittels M Zeitreihen der Laenge L
    const int M = 100;
    const int L = 1000;
    double var  = 0;
    double sum  = 0;
    double varianz[L];
    double zeitreihen[M][L];

    // Die Zeitreihen
    for(int k=0; k<M; k++)
    {
        zeitreihen[k][0] = 1.8;
        zeitreihen[k][1] = 2.0;
        for(int i=2; i<L; i++)
        {
            zeitreihen[k][i] =  a1   * zeitreihen[k][i-1] 
                              + a2   * zeitreihen[k][i-2] 
                              + 0.01 * gasdev(&idum);            
        }
    }

    // Schaetzen der Varianz
    for(int i=0; i<L; i++)
    {
        varianz[i] = 0;
        for(int k=0; k<M; k++)
        {
            varianz[i] += zeitreihen[k][i] * zeitreihen[k][i];
        }
        varianz[i] = varianz[i]/M;
        cout << varianz[i] << endl;
    }

    return 0;
}