/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE FOURIER TRANSFORM OF RIEMANN SPECTRUM			      C
C  06/20/19 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
#include "zero1.h"    // 100000 zeta function zeros
extern char *malloc();
/******************************************************************************
*									      *
*  INTERPOLATE BETWEEN ZEROS						      *
*  04/02/10 (dkc)							      *
*									      *
******************************************************************************/
double interp1(double ND, double ND1, unsigned int I, double init,
	      double *out) {
unsigned int i;
double F,inc;
inc=(ND1-ND)/(double)I;
F=init;
for (i=0; i<I; i++) {
   out[i]=F;
   F=F+inc;
   }
return ND1;
}
//
unsigned int hmax=300; // number of s values
unsigned int MAXN=100000;  // number of zeros used to compute sum
unsigned int I=5;	   // number of values to interpolate between zeros (plus 1)
unsigned int numlim=20000;  // numlim*I must equal MAXN
double factor=0.1; // s value increment, I must exactly divide 1.0/factor
//
void main() {
unsigned int h,i,N;
double *output;
double init,sum,temp,s;
FILE *Outfp;
Outfp = fopen("out2g.dat","w");
output=(double*) malloc(800008);
if (output==NULL) {
   printf("not enough memory \n");
   return;
   }
init=riem[0];
for (i=0; i<(numlim+1); i++)
   init=interp1(riem[i],riem[i+1],I,init,&output[i*I]);
//
for (h=1; h<=hmax; h++) {
   s=factor*(double)h;
   sum=0.0;
   temp=log(s*(double)I);
   for (N=1; N<=MAXN; N++)
     sum=sum+cos(temp*output[N-1]);
   printf("h=%d, s=%llf, sum=%llf \n",h,s,sum);
   if ((s*(double)I)==1.0)
      sum=0.0;
   fprintf(Outfp," %llf \n",-sum);
   }
fclose(Outfp);
return;
}