/* Calculation of Inference Range according a magnetic      */
/* fieldstrength limit, using data from ITU-R PN.368-7      */
/* and from CCIR 322-3.                                     */
/* Version 1.10                                             */

#include <stdio.h>
#include <stddef.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#define pi 3.1415926
#define c 299792458
#define E_asymptote 109.5
#define BW_normal 2700
#define BW_CISPR_lf 200
#define BW_CISPR_hf 9000
#define FAIL 0
#define Success 1
#define TRUE 1
#define FALSE 0
#define INTRO \
"         INTERFERENCE RANGE CALCULATION\n\
                  Version 1.10.\n\n\
              C 2002 T.W.H. Fockens NEDAP.\n\n\
This program calculates the interference range of an \n\
inductive loop system from the measured magnetic field \n\
strength at a given measuring range according to the propagation\n\
model given in the ERC Report 69.\n\n\
For the calculation data about the groundwave propagation\n\
from the ITU-R PN.368-7 has to be inputted as well as the\n\
noise fieldstrength, which has been derived from the ITU-R\n\
Recommendation P.372.\n\
These data is supplied with this program as two graphics. \n\
==========================================================\n\n"

double radiated_power(double freq, double K, double H );
double r,m,K,K_coplanar,K_coaxial, d_transition, d_interference,P_erp_dB,P_erp_nW;
double freq, H_dBuA, H, lab_rad, P_erp, E_1kW, E_int_1km, E_noise, E_noise_BW, MaxINR, Max_E_int, E_int, H_int, BW, BW_ratio;
FILE *fptr;
char filename[13]={"LOGFILE"};
char s1[40],s2[40],s3[40],s4[40];
char answer1[40],answer2[40],answer3[40],answer4[40],answer5[40];
int printfile;
int groundwave;
int broadband;
int NotNoise;
int repeat;


int main(void)
{
  printf(INTRO);
  printf("Do you want to save the result in a logfile, \n");
  printf("which can be imported in a texteditor for printing,(y/n):\n");
  gets(s2);
  if (strlen(s2)) strcpy(answer1,s2);
  if (strcmp(answer1,"y")==0 || strcmp(answer1,"Y")==0) printfile = TRUE;

  if (printfile == TRUE)
  {
    printf("The logfile will be named LOGFILE. \n");
    if ((fptr = fopen(filename, "w"))== NULL)
    {
      fclose( fptr );
      return FAIL;
    }
  }
  if (printfile == TRUE) fprintf(fptr, INTRO);

  do
  {
    printf("Input frequency in kHz: %8.2f\n",freq/1e3);
    gets(s1);
    if (strlen(s1)) freq = atof(s1)*1e3;

	 printf("Input the magnetic field strength limit in dBuA/m: %5.1f\n",H_dBuA);
    gets(s1);
    if (strlen(s1)) H_dBuA = atof(s1);

    printf("Input the measuring distance in meter: %5.1f\n",r);
    gets(s1);
    if (strlen(s1)) r = atof(s1);

    printf("Is the victim receiver at groundlevel (groundwave propagation) (Y),\n");
    printf("or airborne (free space propagation) (N): %s\n",answer3);
    gets(s2);
    if (strlen(s2)) strcpy(answer3,s2);
    if (strcmp(answer3,"y")==0 || strcmp(answer3,"Y")==0) groundwave = TRUE;
    else groundwave = FALSE;

    if (groundwave == TRUE)
    { 
      printf("Input E_1kW@1km according ITU-R PN.368-7: %5.1f\n",E_1kW);
      gets(s1);
      if (strlen(s1)) E_1kW = atof(s1);

      d_transition = 1000*pow(10,(E_1kW - E_asymptote)/20);
      printf("20 to 40 dB/decade transition distance: %4.0f", d_transition ),     printf(" m.\n");
    }
    
    printf("Is the interference range limited by a service defined interference \n level (Y), or by the noise level (N): %s\n",answer5);
    gets(s3);
    if (strlen(s3)) strcpy(answer5,s3);
    if (strcmp(answer5,"y")==0 || strcmp(answer5,"Y")==0) NotNoise = TRUE;
    else NotNoise = FALSE;

    if (NotNoise == TRUE)
    {
		printf("Input the max. acceptable interference level in dBuV/m: %5.1f\n",Max_E_int);
      gets(s1);
      if (strlen(s1)) Max_E_int = atof(s1);
    }
    else
    {
		printf("Input E_noise according ITU-R Rec. P.372 for 2.7 kHz bandwidth in dBuV/m: %5.1f\n",E_noise);
      gets(s1);
      if (strlen(s1)) E_noise = atof(s1);
    }

    printf("Input the bandwidth of the victim receiver in Hz: %8.2f\n",BW);
    gets(s1);
    if (strlen(s1)) BW = atof(s1);
    
    if (freq > 150000)
      BW_ratio = 10*log10(BW_CISPR_hf/BW);
    else
      BW_ratio = 10*log10(BW_CISPR_lf/BW);
      
    printf("Is the interference broadband (Y), or single frequency (N): %s\n",answer4);
    gets(s2);
    if (strlen(s2)) strcpy(answer4,s2);
    if (strcmp(answer4,"y")==0 || strcmp(answer4,"Y")==0) broadband = TRUE;
    else broadband = FALSE;

    if (broadband == TRUE)
    {
      MaxINR = BW_ratio;
      printf("Broadband interference; bandwidth ratio: %5.1f", MaxINR),
      printf("  dB.\n\n");
    }
    else MaxINR = 0;

    if (printfile == TRUE)
    {
      fprintf(fptr, "The input data are: \n");
      fprintf(fptr, "Frequency                              : %8.2f",freq/1e3),
        fprintf(fptr, " kHz.\n");
      fprintf(fptr, "The magnetic field strength limit      : %7.1f", H_dBuA),
		  fprintf(fptr, "  dBuA/m.\n");
      fprintf(fptr, "The measuring distance                 : %7.1f", r),
        fprintf(fptr, "  m.\n");
      if (groundwave == TRUE)
      {
        fprintf(fptr, "The E_1kW@1km according ITU-R PN.368-7 : %7.1f", E_1kW),
			fprintf(fptr, "  dBuV/m.\n") ;
      }
      if (NotNoise == TRUE)
      {
        fprintf(fptr, "The max. acceptable interference level : %7.1f", Max_E_int),
		  fprintf(fptr, "  dBuV/m.\n");
      }
      else
      {
        fprintf(fptr, "The E_noise according ITU-R Rec. P.372 : %7.1f", E_noise),
			 fprintf(fptr, "  dBuV/m.\n");
      }
      fprintf(fptr, "The bandwidth of the victim receiver   : %8.2f", BW),
        fprintf(fptr, " Hz.\n");

      fprintf(fptr, "\n     The results are:\n\n");


      if (groundwave == TRUE)
      {
        fprintf(fptr, "The 20/40 dB/decade transition distance: %5.1f", d_transition),
        fprintf(fptr, "  m.\n"); 
      }
      if (broadband == TRUE) 
        fprintf(fptr, "Broadband interference; bandwidth ratio: %5.1f", MaxINR),
        fprintf(fptr, "  dB.\n\n");
    }
  
    H = 1e-6*pow(10,(H_dBuA/20));
    lab_rad = c/(2*pi*freq);
    K_coaxial = 2*pi*lab_rad*r*r*r/sqrt(lab_rad*lab_rad + r*r);
    K_coplanar = 4*pi*pow(lab_rad,2)*r*r*r/sqrt(pow(lab_rad,4) - lab_rad*lab_rad*r*r + r*r*r*r);
    if (r > lab_rad*2.354)
    {
      K = K_coplanar;
		printf("The field strength at the measuring position is maximal in\n");
      printf("the coplanar direction.\n\n");
      if (printfile == TRUE)
      {
		  fprintf(fptr, "The field strength at the measuring position is maximal in\n");
        fprintf(fptr, "the coplanar direction.\n\n");
      }
    }
    else
    {
      K = K_coaxial;
		printf("The field strength at the measuring position is maximal in\n");
      printf("the coaxial direction.\n\n");
      if (printfile == TRUE)
      {
        fprintf(fptr, "The field strength at the measuring position is maximal in\n");
        fprintf(fptr, "the coaxial direction.\n\n");
      }  
    }

    m = K*H;
    printf("frequency                             :  %8.2f", freq/1e3),
      printf(" kHz.\n");
    printf("Magnetic dipole moment                :      %7.1e", m),
		printf("  A.m2.\n");
    printf("distance                              :  %7.1f", r),
      printf(" m.\n");
    printf("Field strength                        :      %7.1e", H),
      printf(" A/m.\n");
    P_erp = radiated_power(freq, K, H);
    P_erp_dB = 10*log10(P_erp) - 30;
    P_erp_nW = P_erp*1e9;
    printf("Effective radiated power              :  %7.1f", P_erp_dB),
      printf(" dBkW"),
    printf(" / %5.1e", P_erp_nW), printf("  nW.\n");
    
    if (printfile == TRUE)
    {
      fprintf(fptr, "Magnetic dipole momemt                 :     %7.1e",m),
		  fprintf(fptr, "  A.m2.\n");
      fprintf(fptr, "The Effective radiated power           : %7.1f", P_erp_dB),
        fprintf(fptr, "  dBkW"),
        fprintf(fptr, " / %5.1e", P_erp_nW),
        fprintf(fptr, " nW.\n");
    }

    E_int_1km = E_1kW + P_erp_dB;
    if (NotNoise == TRUE)
      E_int = Max_E_int;
    else
    {
      E_noise_BW = E_noise + 10*log10(BW/BW_normal) ;
      E_int = E_noise_BW;
    }
      
    E_int = E_int + MaxINR;
    H_int = pow(10,((E_int - 171.5)/20));
    
    if (groundwave == TRUE)
    {
      d_interference = 1000*pow(10,(E_int_1km - E_int)/40);
/*                   printf("d_interference = %5.0f", d_interference ), printf(" m.\n");           */
      if ((d_interference > d_transition) && (d_interference > lab_rad*2.354)) 
      {
        printf("The interference range extends into the 40 dB/decade range.\n\n");
        printf("The groundwave interference range is  :%5.0f", d_interference ),
          printf(" m.\n");
        if (printfile == TRUE)
        {
          fprintf(fptr, "The interference range extends into the 40 dB/dec. range.\n\n");
          fprintf(fptr, "The groundwave interference range is   : %5.0f", d_interference ),
            fprintf(fptr," m.\n\n");
        }
      }
      else
      {
        d_interference = pow(10,(169.5 + P_erp_dB - E_int)/20);
/*                     printf("d_interference = %5.0f", d_interference ), printf(" m.\n");                */
        if (d_interference > lab_rad*2.354)
        {
          printf("The interference range is limited to the 20 dB/dec. roll-off range.\n\n");
          printf("The groundwave interference range is  : %5.0f", d_interference ),
            printf(" m.\n");
          if (printfile == TRUE)
          {
            fprintf(fptr, "The interference range is limited to the 20 dB/dec. roll-off range.\n\n");
            fprintf(fptr, "The groundwave interference range is   : %5.0f", d_interference ),
              fprintf(fptr," m.\n\n");
          }
        }
        else
        {
          d_interference = sqrt(m/(H_int * lab_rad * 2 * pi));
/*                       printf("d_interference = %5.0f", d_interference ), printf(" m.\n");              */
          if (d_interference > lab_rad)
          {
            printf("The interference range is close to the near field range.\n\n");
            printf("The groundwave interference range is  : %5.0f", d_interference ),
              printf(" m.\n");
            if (printfile == TRUE)
            {
              fprintf(fptr, "The interference range is close to the near field range.\n\n");
              fprintf(fptr, "The groundwave interference range is   : %5.0f", d_interference ),
                fprintf(fptr," m.\n\n");
            }
          }
          else
          {
            d_interference = pow((m/(2 * pi * H_int)),(1.0/3));
/*                         printf("d_interference = %5.0f", d_interference ), printf(" m.\n");              */
            printf("The interference range is inside the near field range.\n\n");
            printf("The groundwave interference range is   : %5.0f", d_interference ),
              printf(" m.\n");
            if (printfile == TRUE)
            {
              fprintf(fptr, "The interference range is inside the near field range.\n\n");
              fprintf(fptr, "The groundwave interference range is  : %5.0f", d_interference ),
                fprintf(fptr," m.\n\n");
            }
          }

        }
      }
    }
    else
    {
      d_interference = pow(10,(169.5 + P_erp_dB - E_int)/20);
/*                     printf("d_interference = %5.0f", d_interference ), printf(" m.\n");                */
      if (d_interference > lab_rad*2.354)
      {
        printf("The free space interference range is  : %5.0f", d_interference ),
          printf(" m.\n");
        if (printfile == TRUE)
        {
          fprintf(fptr, "The free space interference range is   : %5.0f", d_interference ),
            fprintf(fptr," m.\n\n");
        }
      }
      else
      {
        d_interference = sqrt(m/(H_int * lab_rad * 2 * pi));
/*                       printf("d_interference = %5.0f", d_interference ), printf(" m.\n");              */
        if (d_interference > lab_rad)
        {
          printf("The interference range is close to the near field range.\n\n");
          printf("The free space interference range is  : %5.0f", d_interference ),
            printf(" m.\n");
          if (printfile == TRUE)
          {
            fprintf(fptr, "The interference range is close to the near field range.\n\n");
            fprintf(fptr, "The free space interference range is   : %5.0f", d_interference ),
              fprintf(fptr," m.\n\n");
          }
        }
        else
        {
          d_interference = pow((m/(2 * pi * H_int)),(1.0/3));
/*                         printf("d_interference = %5.0f", d_interference ), printf(" m.\n");              */
          printf("The interference range is inside the near field range.\n\n");
          printf("The free space interference range is  : %5.0f", d_interference ),
            printf(" m.\n");
          if (printfile == TRUE)
          {
            fprintf(fptr, "The interference range is inside the near field range.\n\n");
            fprintf(fptr, "The free space interference range is   : %5.0f", d_interference ),
              fprintf(fptr," m.\n\n");
          }
        }
      }
    }

    printf("=================================================================\n");
      if (printfile == TRUE)  
        fprintf(fptr, "=======================================================================\n\n");
    printf("Do you want to do another calculation (Y), or to exit (N): %s\n",answer2);
    gets(s2);
    if (strlen(s2)) strcpy(answer2,s2);
    if (strcmp(answer2,"y")==0 || strcmp(answer2,"Y")==0) repeat = TRUE;
    else repeat = FALSE;
  }
  while  (repeat==TRUE);
  fclose(fptr);
  return Success;
}

double radiated_power(double freq, double K, double H )
{
        double  Result;
        Result = 3.848e-30*pow(freq,4)*K*K*H*H;
        return(Result);     
}