#ifndef __RANDNUM_H99_4_5_22080409
#define __RANDNUM_H99_4_5_22080409

#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */

float rand1(void); /* Uniform deviation within [0, 1] */
float ran1(long *idum); /* Uniform deviation within [0, 1] */

#ifdef __cplusplus
}
#endif /* __cplusplus */

#endif /* __RANDNUM_H99_4_5_22080409 */

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

#include <math.h>
#include <time.h>
#include <stdlib.h> 
#include <iostream.h>
#include <fstream.h>

static long random_number_seed;
int flag_seeded=0;

float rand1(void)  /* Generating Uniform random number between 0.0 and 1.0*/
{
        if(!flag_seeded)
        {
                random_number_seed = (long)time(NULL);
                random_number_seed = 1664525L*random_number_seed + 1013904223L;
                if(random_number_seed>0)
                        random_number_seed = -random_number_seed;
                flag_seeded=1;
        }
        return ran1(&random_number_seed);
}

float ran1(long *idum)
{
        int j;
        long k;
        static long iy=0;
        static long iv[NTAB];
        float temp;
        if (*idum <= 0 || !iy)
        {                                                        /*Initialize.*/
                if (-(*idum) < 1)
                        *idum=1;                                 /*Be sure to prevent idum = 0.*/
                else
                        *idum = -(*idum);
                for (j=NTAB+7;j>=0;j--)
                {                                                /*Load the shuffle table (after 8 warm-ups).*/
                        k=(*idum)/IQ;
                        *idum=IA*(*idum-k*IQ)-IR*k;
                        if (*idum < 0)
                                *idum += IM;
                        if (j < NTAB)
                        iv[j] = *idum;
                }
                iy=iv[0];
        }
        k=(*idum)/IQ; /* Start here when not initializing. */
        *idum=IA*(*idum-k*IQ)-IR*k; /* Compute idum=(IA*idum) % IM without overflows */
        if (*idum < 0)
                *idum += IM;            /* by Schrage's method. */
        j=iy/NDIV; /* Will be in the range 0..NTAB-1. */
        iy=iv[j];  /* Output previously stored value and rell the shuffle table. */
        iv[j] = *idum;
        if ((temp=AM*iy) > RNMX)
                return RNMX; /* Because users don't expect endpoint values. */
        else
                return temp;
}

void main()
{
     ofstream of("dtate.txt");         
     fstream iosp;
	 iosp.open("spectral.dat",ios::out|ios::binary|ios::in);
     float d,w,dw;
     float *res;
     res=new float [201*100];
     dw=100.0/1000.0;
     float t;
     float fi;
     float s=0.0005;      
     for(int i=0;i<201*100;i++)
     {    
          if(i%500==0)
               cout<<i<<endl;
          res[i]=0.0;
          t=i*0.0314159265358979/2.0;
          for(int j=0;j<1000;j++) ///////
          {
              //int rnd=rand()%2500;
              //d=(float)(rnd)/2500.0;
              //d=0.0005*d-0.0005/2.0;
              w=-50.0+(j-0.5)*dw;
              fi=2.0*3.1415926*rand1();
              res[i]=res[i]+sqrt(2.0*s*dw)*cos(w*t+fi);
          }
         /// of<<t<<"             "<<res[i]<<endl;		  
     }
	 float temp;
     for(i=0;i<201*100;i++)
	 {
		 temp=res[i];
		 iosp.write((char *)(&temp),sizeof(float));
	 }
	 float rt;
	 for(i=0;i<201;i++)
	 {
		 rt=0.0;
		 for(int j=0;j<8500;j++)
		 {
			 rt=rt+res[0+j]*res[i+j];
		 }
		 rt=rt/8500.0;
		 of<<i*0.0314159265358979/2.0<<"      "<<rt<<endl;

	 }

	 
     //////////////////////////////////////////////////////     
     /*
	 float rt[100];     
     for(i=0;i<100;i++)
     {
          cout<<i<<endl;
          t=i*0.125;
          rt[i]=0.0;
		  int num=0;
		  float temp;
		  temp=0.0;
          for(int k=0;k<800;k++)
          {
             for(int j=k;j<800;j++)
             {
                int te=j-k;
                if(te<0)
                { 
                     te=-te;
                }
                if(te==i)/////////////////////time ///////////k 與 j冽相乘 
                {
					  num++;
					  
					  //temp=0.0;
                      for(int r=0;r<250;r++)
                      {
                            temp=temp+res[k+r*800]*res[j+r*800];
                      }
                      //temp=temp/250.0;
					  //rt[i]=(rt[i]+temp)/2.0;
                }
             }
          }
		  rt[i]=temp/(num*250);
          of<<t<<"         "<<rt[i]<<endl;
     }*/
       delete [] res;
	   
    /*
     float real[500];
     float ima[500];  
     for(int i=0;i<500;i++)
     {
         w=-50+i*0.2;
         real[i]=0;
         ima[i]=0;
         for(int j=0;j<256;j++)
         {
             t=0.005*j;
             real[i]=real[i]+rt[j]*cos(w*t)*0.005;
             ima[i]=ima[i]+rt[j]*sin(w*t)*0.005;
         }
         real[i]=real[i]/(2*3.1415926);
         ima[i]=ima[i]/(2*3.1415926);     
     } 
     
     
     float tre[500];
     for(int i=0;i<500;i++)
     {
         tre[i]=sqrt(real[i]*real[i]+ima[i]*ima[i]);
         w=-50+i*0.2;
         of<<w<<"             "<<tre[i]<<endl;
     }
     */      
     
   
               
}