#include<stdio.h>
#include<math.h>
#include<conio.h>

//variables introduced by the user
int Sprime;                      
int Lprime;
int k;
int N;
double minint;

// variables recycled in various parts of the algorithm
int j;                       
int minimal;
double DeltaD;
char yn;

//varible used to number the reflections in the basic structure
double t;

//value Nb for the basic structure
double Nb;

// array containing the intensity values. assumes that N<=1000 
double intensity[1000];      

//variable used to number the reflections in the derived structure
int nr;                     

//numbers used to calculate if S,L coprime
int a,b,r,q;    

//variables used to determine when a basic structure with Nb respects the minimal shift condition
int ding, refl;       

//definition for pi number
#define PI 3.14159265358979323846            

//define file variable 
FILE  *fv;

// Function cmd(S',L') calculates if S', L' are coprime numbers using an Euclidean algorithm and turns them into coprime numbers if they are not
void cmd(int *Sprime, int *Lprime)
{
  if(*Sprime>*Lprime)
    {
      a=*Sprime;
      b=*Lprime;
    }
  else
    {
      a=*Lprime;
      b=*Sprime;
    }
  do
    {
      q=a/b;
      r=a-b*q;
      if(r==0)
	break;
      else
	{
	  a=b;
	  b=r;
	}
    }while(r!=0);
  *Sprime=*Sprime/b;
  *Lprime=*Lprime/b;
}


// Function I(D) calculates the intensity due to shift D using the (sin(Pi*D)/(Pi*D))^2 formula
double I(double D)
{ 
  //to avoid division by 0 if D is very small the program returns 1 which is the analytical value for the function
  if(D<1E-50)
    return 1;
  else
    return (sin(PI*D)/(PI*D))*(sin(PI*D)/(PI*D));
}

int main(void)
{ 
  fv=fopen("Intensities.txt", "w");
  
  beginning: 	
  
  //Prompts the user to enter the values S',L' and k
  do
    {
      printf("Please enter the values for S', L' and k. \n");
      printf("S'=");
      fflush(stdout);
      scanf("%d", &Sprime);
      printf("L'=");
      fflush(stdout);
      scanf("%d", &Lprime);
      printf("k=");
      fflush(stdout);
      scanf("%d", &k);
      if((Sprime<0)||(Lprime<0))
	printf("S' and L' must be positive. \n");
    }while ((Sprime<0)||(Lprime<0));
  
  //Sequence prompting the introduction of a threshold value for the intensities to be displayed 
  printf("Please enter the minimal threshold value for the intensities to be displayed. Default value is 0.09 \n");
  printf("Minimal intensity=");
  fflush(stdout);
  scanf("%lf", &minint); 
  
  //Tests whether S' or L' are 0. If either of them is 0 then the non-zero value is converted into 1 since this sequence is equivalent to a sequence with that value equal to 1
  if((Sprime!=0)&&(Lprime!=0))         
    cmd(&Sprime,&Lprime);              
  else                      
    {
      if(Sprime==0)
      Lprime=1;
      else Sprime=1;
    }
  
  //Calculate and display N
  N=Sprime*(2*k+3)+Lprime*(2*k+5);      
  printf("Calculated N= %d \n",N);
  
  //Writes the values S', L', k and N to the file Intensities.txt
  fprintf(fv, "S'= %d \n", Sprime);
  fprintf(fv, "L'= %d \n", Lprime); 
  fprintf(fv, "k = %d \n", k);
  fprintf(fv, "N = %d \n", N);
  
  //Calculates position 'minimal' in the pattern where the minimal shift condition applies
  minimal=(N-Sprime-Lprime)/2;          
                              
  //Calculates a starting Nb value for the basic structure
  Nb=2*k+3;                   
  
  //main body of the program which searches for the basic structure Nb and calculates the intensity array for each Nb
  do
    { 
      refl=0;
      ding=0;      
      
      //Resets all the values in the intensity array before each Nb is considered
      for(nr=0;nr<=N;nr++)     
	{
	  intensity[nr]=0;
	}  
           
      //Sequence calculating the intensity array for a particular Nb and testing if the minimal shift condition is satisfied
      for(nr=0; nr<=N; nr++)                   
	{ 
	  for(t=0; t<=Nb; t++)
	    {
	      //Calculates the shift between position nr and position t
	      DeltaD=fabs((t/Nb)*N-nr);      
	      		
	      
	      //The value of DeltaD-1/Nb is compared with an arbirarily small number 1E-10 to determine if they are equal
	      if (fabs((DeltaD)-(1/Nb))<=1E-10)  
		{ 
		  //If the condition if satisfied for other reflections than those satisfying the minimal shif condition then the program changes the variable refl and discards the base
		  if ((nr!=minimal)&&(nr!=minimal+Sprime+Lprime))
		    {refl++;}
		  
		  //If the minimal shift condition is satisfied the variable ding is used to specify this base is correct
		  if ((nr==minimal)&&(Nb!=N))
		    {ding=1;}                      
		}                              
		
	      //Sums the intensities for a particular nr
	      intensity[nr]=intensity[nr]+I(DeltaD);
	    }
	}
      //Verifies if the base Nb satisfies the minimal shift condition by checking the variable ding and refl
      if((ding==1)&&(refl==0))                              
	{                                      
	  //intensities greater than the value of 'minint' are displayed along with their respective nr
	  fprintf(fv,"\n Base Nb=  %.0lf \n", Nb);
	  printf("\n Base Nb=  %.0lf \n", Nb);        
	  fprintf(fv,"nr     Intensity \n");
	  printf("nr     Intensity \n");
	  for(j=0;j<=N;j++)
	    {
	      if(intensity[j]>minint)
		{
		   fprintf(fv,"%d %12.5lf\n", j, intensity[j]);
	           printf("%d %12.5lf\n", j, intensity[j]);    
		}                                            
	    }
	}
      
      //Increments Nb by 1
      Nb++;   
    }while(Nb<N+3);
  
  fprintf(fv, "\n \n");
  
  //Sequence which checkes if the user wishes to run further simulations in this instance of the program 
  do
  {
     printf("Do you wish to run another simulation? (y or n) : ");
     fflush(stdout);
     scanf("%s", &yn);
     if (yn=='y')
     	{goto beginning;}	
     else {if (yn=='n') goto end;}
     fclose(fv);
  }while ((yn!='n')&&(yn!='y'));
  
  end:
  printf("\nThe intensities diplayed on the screen are also written in the text file Intensities.txt located in the current folder where the program is being run. \n");
  printf("Please press enter to quit.");
  getch();
  return 0;
}