//SIDI11II.C   imaging with SIDI 1.1   takoj narise sliko

//             Marko Cebokli  mar 2005   GNU/GPL licence

//             www.gnu.org/licences/licences.html#GPL

//             software sample timing

//input connection LPT port D25
//				   pin 11 channel 1 I             bit 7
//                 pin 10 channel 1 Q             bit 6
//                 pin 12 channel 2 I             bit 5
//                 pin 13 channel 2 Q             bit 4
//                 pin 18..25  ground

//daj se:

//  - ce synt ne locka itd  (abort)
//  - C: D: itd
//  - vprasaj za avg

#include <stdio.h>
#include <time.h>
#include <dos.h>       //inportb
#include <conio.h>     //kbhit,getch
#include <stdlib.h>    //exit()
#include <math.h>
#include <graphics.h>
#include <string.h>
#include <dir.h>

void cddd(void);
void cddg(void);
void cgdd(void);
void cgdg(void);
void i2creset(void);
int  prog3302(unsigned char bb[]);
int  beri1935(unsigned char adr);
int  sl1935(int adr, double frek, double *pfrek, int *rdiv, int *mdiv);
void mkcortab(char faza,char real[],char imag[]);
void mkdctab(char dctab[]);
void vzorcenje(void);
void checkdc(void);
void checkquad(void);
void korelator(float *korr,float *kori);
void checkspeed(float *brate, float *srate);
void nastavi_uro(void);
void pitaj_budalu(char sname[], char comment[], char assfile[]);
void risi_franze(float brate, float srate, double pfrek, int rdiv1, int mdiv1, int rdiv2, int mdiv2);
void risi_polar();
void pisi_korel();
int  open_datax();
void odrezi(char *string);
void pisi_meta(double ts, float brate, float srate, double pfrek, int rdiv1, int mdiv1, int rdiv2, int mdiv2, char sname[], char comment[], char assfile[]);
void set_freq(double *frek, double *pfrek, int *rdiv1, int *mdiv1, int *rdiv2, int *mdiv2);
void tests(float *brate, float *srate);
void fft(unsigned int n, float xr[], float xi[], float wr[], float wi[], int bi[]);
void fftmojprip(unsigned int n,float wr[], float wi[], int bi[]);
void setpal();
void risisl(int vm, float din);
void fft_2D(int vm);
void utez(int vm, int tip);
void promuckaj(int vm);
void cakitipko();


int LPT=0x378;
float PI=3.141592654;

unsigned int sig[32767];
char ct0real[65535],ct0imag[65535]; //veliki arrayi zunaj, sicer crash
char dctab[65535];

float wr[33],wi[33];
int bi[33];
float xr[33],xi[33];
float slikar[32][32],slikai[32][32];

FILE *ff,*fm;       //bile so tezave, ce sem jih prenasal kot parameter

char irqmask;

//----------------------------------------------------------------------
//------------routines for tuner synth (SL1935) programming-------------

//d0 (pin2) = SCL,  d1 (pin3) = SDA,  error (pin15) = vhod SDA

void cddd()                        //SCL low,  SDA low  (transistor invert!)
{outportb(LPT,0xFF);outportb(LPT,0xFF);outportb(LPT,0xFF);outportb(LPT,0xFF);
outportb(LPT,0xFF);outportb(LPT,0xFF);outportb(LPT,0xFF);outportb(LPT,0xFF);
outportb(LPT,0xFF);outportb(LPT,0xFF);outportb(LPT,0xFF);outportb(LPT,0xFF);}

void cddg()                        //SCL low,  SDA high
{outportb(LPT,0xFD);outportb(LPT,0xFD);outportb(LPT,0xFD);outportb(LPT,0xFD);
outportb(LPT,0xFD);outportb(LPT,0xFD);outportb(LPT,0xFD);outportb(LPT,0xFD);
outportb(LPT,0xFD);outportb(LPT,0xFD);outportb(LPT,0xFD);outportb(LPT,0xFD);}

void cgdd()                        //SCL high, SDA low
{outportb(LPT,0xFE);outportb(LPT,0xFE);outportb(LPT,0xFE);outportb(LPT,0xFE);
outportb(LPT,0xFE);outportb(LPT,0xFE);outportb(LPT,0xFE);outportb(LPT,0xFE);
outportb(LPT,0xFE);outportb(LPT,0xFE);outportb(LPT,0xFE);outportb(LPT,0xFE);}

void cgdg()                        //SCL high, SDA high
{outportb(LPT,0xFC);outportb(LPT,0xFC);outportb(LPT,0xFC);outportb(LPT,0xFC);
outportb(LPT,0xFC);outportb(LPT,0xFC);outportb(LPT,0xFC);outportb(LPT,0xFC);
outportb(LPT,0xFC);outportb(LPT,0xFC);outportb(LPT,0xFC);outportb(LPT,0xFC);}

void i2creset(void)
{cddg();cgdg();cgdd();cgdg();}

//--------------------------------------------------------------------

int prog3302(unsigned char bb[])        //send 5 bytes bb[] to I2C bus
{
int i,j,a;unsigned char b;
cgdd();cddd();   //START
a=0;
for (i=0;i<5;i++)
  {
  b=bb[i];
  for (j=0;j<8;j++)
	{
	if ((b&128)==128)
	  {cddg();cgdg();cddg();}
	else
	  {cddd();cgdd();cddd();}
	b=b<<1;
	}
  cddg();cgdg();a=a+(inportb(LPT+1)&8);cddg();	  //ACK
  }
cgdd();cgdg();        //STOP
return a/8;
}

//-----------------------------------------------------------------

int beri1935(unsigned char adr)       //read status byte from sl1935
{
int i,a;unsigned char b;
cgdd();cddd();        //START
b=192+2*adr+1;
for (i=0;i<8;i++)
  {
  if ((b&128)==128)
	{cddg();cgdg();cddg();}
  else
	{cddd();cgdd();cddd();}
  b=b<<1;
  }
cddg();cgdg();a=(inportb(LPT+1)&8);cddg();	  //ACK
b=0;
for (i=0;i<8;i++)
  {
  cddg();cgdg();b=(b<<1)+((inportb(LPT+1)&8)/8);cddg();
  }
cddg();cgdg();cddg();     //ACK
cddd();cgdd();cgdg();     //STOP
if (a==0) return b; else return -a;
}

//---------------------------------------------------------------
//function sl1935 programs the Zarlink SL1935 synthesizer
//adr = address 0..3,  frek MHz
//if successful returns 0

int sl1935(int adr, double frek, double *pfrek, int *rdiv, int *mdiv)
{
unsigned char bb[5],cp,bp,buf,tm,rr,vco,p0;
double xtal,fref;
unsigned int n;
int stat,statb;
int rd[]={2,4,8,16,32,64,128,256,0,5,10,20,40,80,160,320,0,
6,12,24,48,96,192,384,0,7,14,28,56,112,224,448}; //ref divider ratios

cp=0;   //charge pump  0...3  (155,330,690,1450 uA)
bp=1;   //baseband path 0...1
buf=0;  //BUFREF output 0...1   (1=ON)
tm=0;   //test mode 0...7
rr=3;   //ref divider 0...31 (see rd array)
vco=0;  //VCO select 0...1  (automatic switch at 1390MHz)
p0=0;   //p0 port output 0...1

bb[0]=192+2*adr;
bb[1]=0x3;bb[2]=0xE8;
bb[3]=128+cp*32+rr;
bb[4]=tm*32+vco*16+bp*8+buf*2+p0;

xtal=4000000.00;
xtal=8998500.00;
fref=xtal/rd[rr];

n=frek*1000000.0/fref;
bb[1]=n/256;bb[2]=n%256;
if (frek<1390) bb[4]=bb[4]|16; else bb[4]=bb[4]&239;
printf("\n @%1d ",adr);
printf(" F= %5.6f MHz ",n*fref/1000000.0);
printf("(Fref= %2.6f, N= %6d )",fref/1000000.0,n);
stat=prog3302(bb);
if (stat!=0) printf("\n Prog3302: no ACK, stat= %d",stat);
sleep(1);
statb=beri1935(adr);
if (statb<0)
  printf("\n Beri1935: no ACK");
else
  {printf(" POR: %d",(statb&128)/128);printf("  LOCK: %d",(statb&64)/64);}

*pfrek=n*fref/1000000.0;
*rdiv=rd[rr];
*mdiv=n;
return stat+abs(statb);
}


//--------------------------------------------------------
//mkcortab prepares lookup tables for the correlator
//data input format: AaBbAaBbAaBbAaBb  (chans A and B,  In phase X quadrature x)
//can prepare the table for 4 possible phases - faza=0,1,2,3
//0=0deg, 1=90deg, 2=180 deg, 3=270deg (for fringe rotation)

void mkcortab(char faza,char real[],char imag[])
{
char faza0r[]={1,0,0,-1,0,1,-1,0,0,-1,1,0,-1,0,0,1};
char faza0i[]={0,1,-1,0,-1,0,0,1,1,0,0,-1,0,-1,1,0};
char faza90r[]={0,1,-1,0,-1,0,0,1,1,0,0,-1,0,-1,1,0};
char faza90i[]={-1,0,0,1,0,-1,1,0,0,1,-1,0,1,0,0,-1};
char faza180r[]={-1,0,0,1,0,-1,1,0,0,1,-1,0,1,0,0,-1};
char faza180i[]={0,-1,1,0,1,0,0,-1,-1,0,0,1,0,1,-1,0};
char faza270r[]={0,-1,1,0,1,0,0,-1,-1,0,0,1,0,1,-1,0};
char faza270i[]={1,0,0,-1,0,1,-1,0,0,-1,1,0,-1,0,0,1};
long j;
int b1,b2,b3,b4;

for (j=0;j<65536;j++)
  {
  b1=(j&0x0000000F);
  b2=(j&0x000000F0)>>4;
  b3=(j&0x00000F00)>>8;
  b4=(j&0x0000F000)>>12;
  switch (faza)
	{
	case 0:{real[j]=faza0r[b1]+faza0r[b2]+faza0r[b3]+faza0r[b4];
		   imag[j]=faza0i[b1]+faza0i[b2]+faza0i[b3]+faza0i[b4];break;}
	case 1:{real[j]=faza90r[b1]+faza90r[b2]+faza90r[b3]+faza90r[b4];
		   imag[j]=faza90i[b1]+faza90i[b2]+faza90i[b3]+faza90i[b4];break;}
	case 2:{real[j]=faza180r[b1]+faza180r[b2]+faza180r[b3]+faza180r[b4];
		   imag[j]=faza180i[b1]+faza180i[b2]+faza180i[b3]+faza180i[b4];break;}
	case 3:{real[j]=faza270r[b1]+faza270r[b2]+faza270r[b3]+faza270r[b4];
		   imag[j]=faza270i[b1]+faza270i[b2]+faza270i[b3]+faza270i[b4];break;}
	default:{printf("\n Napaka v mkcortab: faza= %d",faza);exit(1);}
	}
  }
}

//-----------------------------------------------------------
//mkdctab prepares a lookup table for fast DC ofset determination
//tests bits 0,4,8,12   others are done by shifting

void mkdctab(char dctab[])
{
long j;
for (j=0;j<65536;j++)
  {
  dctab[j]=-2;
  if ((j&0x1)==0x1) dctab[j]=dctab[j]+1;
  if ((j&0x10)==0x10) dctab[j]=dctab[j]+1;
  if ((j&0x100)==0x100) dctab[j]=dctab[j]+1;
  if ((j&0x1000)==0x1000) dctab[j]=dctab[j]+1;
  }
}
//-----------------------------------------------------------

void vzorcenje()     //this one does the sampling
{                    //array is written each time to equalize timing
int i;
long j;
unsigned char a;

for (j=0;j<131072;j++)  //131072
  {
  i=(int)(j>>2);
  a=inportb(LPT+1)>>4;    //D25:   *pin11  pin10  pin12  pin13  (pin15)
  sig[i]=(sig[i]<<4)+a;
  }
}

//-----------------------------------------------------------

void checkdc()      //calculate the DC offsets for adjustment
{
long j,dc1,dc2,dc3,dc4;
float dc1p,dc2p,dc3p,dc4p;
unsigned int i;

dc1=0;dc2=0;dc3=0;dc4=0;
for (j=0;j<32768;j++)
  {
  i=sig[j]^0x8888;
  if ((i&0x1)==0x1) dc1=dc1+1;     //Q2
  if ((i&0x2)==0x2) dc2=dc2+1;     //I2
  if ((i&0x4)==0x4) dc3=dc3+1;     //Q1
  if ((i&0x8)==0x8) dc4=dc4+1;     //I1
  if ((i&0x10)==0x10) dc1=dc1+1;
  if ((i&0x20)==0x20) dc2=dc2+1;
  if ((i&0x40)==0x40) dc3=dc3+1;
  if ((i&0x80)==0x80) dc4=dc4+1;
  if ((i&0x100)==0x100) dc1=dc1+1;
  if ((i&0x200)==0x200) dc2=dc2+1;
  if ((i&0x400)==0x400) dc3=dc3+1;
  if ((i&0x800)==0x800) dc4=dc4+1;
  if ((i&0x1000)==0x1000) dc1=dc1+1;
  if ((i&0x2000)==0x2000) dc2=dc2+1;
  if ((i&0x4000)==0x4000) dc3=dc3+1;
  if ((i&0x8000)==0x8000) dc4=dc4+1;
  }
dc1p=(dc1-65536)/655.36;dc2p=(dc2-65536)/655.36;
dc3p=(dc3-65536)/655.36;dc4p=(dc4-65536)/655.36;


//printf("\nDC:  I1=%lu  ",dc4);printf("%lu  Q1=",dc3);
//printf("%lu  I2=",dc3);printf("%lu  Q2=",dc1);
//printf("\n DC offsets: ");
//printf(" I1=%5.2f%%  ",dc4p);printf("Q1=%5.2f%%  ",dc3p);
//printf("I2=%5.2f%%  ",dc2p);printf("Q2=%5.2f%%  ",dc1p);

printf(" %+4.2f%%  %+4.2f%% ",dc4p,dc3p);
printf(" %+4.2f%%  %+4.2f%% ",dc2p,dc1p);

}

//-----------------------------------------------------------

void checkquad()      //check quadrature of I/Q channels
{                               //with noise input, perfect quadrature
long j,dc1,dc2,dc3,dc4,k1,k2;   //means no correlation between I and Q
unsigned int i,k;
float q1,q2;

k1=0;k2=0;dc1=0;dc2=0;dc3=0;dc4=0;
for (j=0;j<32768;j++)
  {
  i=sig[j]^0x8888;k=i;
  k=k^(k>>1);
  if((k&0x1)==0x1) k1=k1+1;
  if((k&0x4)==0x4) k2=k2+1;
  if((k&0x10)==0x10) k1=k1+1;
  if((k&0x40)==0x40) k2=k2+1;
  if((k&0x100)==0x100) k1=k1+1;
  if((k&0x400)==0x400) k2=k2+1;
  if((k&0x1000)==0x1000) k1=k1+1;
  if((k&0x4000)==0x4000) k2=k2+1;
  dc1=dc1+(long)dctab[i];            //Q2
  i=i>>1;dc2=dc2+(long)dctab[i];     //I2
  i=i>>1;dc3=dc3+(long)dctab[i];     //Q1
  i=i>>1;dc4=dc4+(long)dctab[i];     //I1
  }
q1=(k1-65536.0)/65536;   //tukaj se DC korekcija!
q2=(k2-65536.0)/65536;

//printf("\n  Quadrature %1.5f %1.5f",q1,q2);printf("   %ld  %ld",k1,k2);
printf("   %+2.2f  %+2.2f",180*asin(q1)/PI,180*asin(q2)/PI);
}

//-----------------------------------------------------------

void korelator(float *korr,float *kori)		//here correlation is done using lookup tables
{
long j,rkor,ikor,dc1,dc2,dc3,dc4;
float dckorr,dckori,kr,ki;
unsigned int i;

rkor=0;ikor=0;dc1=0;dc2=0;dc3=0;dc4=0;
for (j=0;j<32768;j++)
  {
  i=sig[j]^0x8888;      //correction for LPT inverted status bits
  rkor=rkor+(long)ct0real[i];
  ikor=ikor+(long)ct0imag[i];
  dc1=dc1+(long)dctab[i];            //Q2
  i=i>>1;dc2=dc2+(long)dctab[i];     //I2
  i=i>>1;dc3=dc3+(long)dctab[i];     //Q1
  i=i>>1;dc4=dc4+(long)dctab[i];     //I1
  }
kr=(float)rkor/131072;
ki=(float)ikor/131072;

dckorr=((float)dc1*dc3+(float)dc2*dc4)/8589934592.0;  //false correlation
dckori=((float)dc3*dc2-(float)dc1*dc4)/8589934592.0;  //caused by DC bias
														//858....  = 2^33
kr=kr-dckorr;
ki=ki-dckori;

kr=sin(PI*kr/2);ki=sin(PI*ki/2);     // VAN VLECK

*korr=kr;*kori=ki;
}

//----------------------------------------------------------

void checkspeed(float *brate, float *srate) //measure speed of computer (sample rate etc)
{
//clock_t sta,mid,sto;
int sta,mid,sto;
float korr,kori;
int i,j;

printf("\n\n\n Measuring speed of computer, please wait cca 15 seconds \n");
//printf("\n clktck=%8.2f \n",CLK_TCK);
i=1;
do
  {
  sta=(int)clock();
  for (j=1;j<i;j++) vzorcenje();            //do the sampling
  mid=(int)clock();
  for (j=1;j<i;j++) korelator(&korr,&kori); //do the correlation
  sto=(int)clock();
//  printf("\n I,tv,tc= %d",i);printf(" %d %d",mid-sta,sto-sta);
  printf("\n I,t= %d %d",i,sto-sta);//printf(" %d %d %d",sta,mid,sto);
  i=2*i;
  }
while ((sto-sta)<100);
i=i/2;

printf("\n Iteracij= %d",i);
printf("\n cas vzorcenja= %d",mid-sta);
printf("\n cas korelacije = %d",sto-mid);
*srate=131072.0/(mid-sta)*CLK_TCK*i;
printf("\n sampling rate= %8.0f /sec",*srate);
*brate=1.0/(sto-sta)*CLK_TCK*i;
printf("\n block rate= %3.2f /sec \n",*brate);

}

//----------------------------------------------------------------

 void nastavi_uro(void)     // set time/date from HW clock
 {
 union REGS regs;

 regs.h.ah=2;                   //BIOS: get current time
 int86(0x1A, &regs, &regs);

 regs.h.ch=(regs.h.ch&0xF)+10*(regs.h.ch>>4);   //BCD > BIN
 regs.h.cl=(regs.h.cl&0xF)+10*(regs.h.cl>>4);
 regs.h.dh=(regs.h.dh&0xF)+10*(regs.h.dh>>4);
 regs.h.dl=0;

 regs.h.ah=0x2D;				//DOS: set system time
 int86(0x21, &regs, &regs);       //Nastavi tudi HW uro!!!

 regs.h.ah=4;                   //BIOS: get current date
 int86(0x1A, &regs, &regs);

 regs.h.ch=(regs.h.ch&0xF)+10*(regs.h.ch>>4);   //BCD > BIN
 regs.h.cl=(regs.h.cl&0xF)+10*(regs.h.cl>>4);
 regs.x.cx=100*regs.h.ch+regs.h.cl;
 regs.h.dh=(regs.h.dh&0xF)+10*(regs.h.dh>>4);
 regs.h.dl=(regs.h.dl&0xF)+10*(regs.h.dl>>4);

 regs.h.ah=0x2B;				//DOS: set system date
 int86(0x21, &regs, &regs);

}

//-----------------------------------------------------------
//ask the user for various data

void pitaj_budalu(char sname[], char comment[], char assfile[])
{
char filn[50],kjesmo[50],string[80];
char *naslov="SIDI11 METADATA";
int imaname;

do
  {
  printf("\n Enter name of file with meta data: ");
  gets(filn);imaname=sscanf(filn,"%30s",filn);fflush(stdin);
  if (!(imaname==1))
	{
	printf("\n no filename, exiting");
	outportb(0x21,irqmask);      //resume timer interrupt
	nastavi_uro();
	exit(1);
	}
  fm=fopen(filn,"r");
  if (fm==NULL)
	{
	printf("\n cannot open file %s\\%s\n",kjesmo,filn);
	}
  else
	{
	fgets(string,80,fm);    //top line, must be "SIDI11 METADATA"
	if (!(memcmp(string,naslov,15)==0))
	  {
	  printf("\n %s\\%s is not a valid SIDI v1.1 metadata file\n",kjesmo,filn);
	  fclose(fm);fm=NULL;
	  }
	}
  }
while (fm==NULL);

printf("\n Enter name of observed source: ");
gets(sname);fflush(stdin);
printf("\n Enter comment: ");
gets(comment);fflush(stdin);
printf("\n Enter name(s) of associated files: ");
gets(assfile);fflush(stdin);

}

//------------------------------------------------------------
//function risi_franze draws fringes and stores correlation values
//with timestamps into a DataX output file

void risi_franze(float brate, float srate, double pfrek, int rdiv1, int mdiv1, int rdiv2, int mdiv2)
{
int i,x,y,avg,size;
float korr,kori,kavgr,kavgi;
time_t t;
struct tm *utc;
double ts;
char sname[80],comment[80],assfile[80];
char string[80],yes='y',*seeno;

float yyy=2134.342;int uuu;

avg=10;                   	//averaging

pitaj_budalu(sname,comment,assfile);

printf("\n Daj velikost matrike ");scanf("%d",&size);fflush(stdin);

do
  {

nastavi_uro();
t = time(NULL);
utc = localtime(&t);
ts=(double)t;

open_datax();
pisi_meta(ts, brate/avg, srate, pfrek, rdiv1, mdiv1, rdiv2, mdiv2, sname, comment, assfile);
fprintf(ff,"\n%d square imaging matrix",size);

clrscr();

printf("                          SIDI 1.1.1 IMAGER by s57uuu  \n\n");

for (x=0;x<size;x++)
  {
  for (y=0;y<size;y++)
	{
	printf("\n Daj v vrsto %d  polozaj %d in pritisni tipko!",x,y);
//	getch();
	cakitipko();
	sound(500);delay(200);nosound();
	nastavi_uro();
	t = time(NULL);
	utc = localtime(&t);
	ts=(double)t;

	kavgr=0.0;kavgi=0.0;
	for (i=0;i<avg;i++)
	  {
	  vzorcenje();
	  korelator(&korr,&kori);
	  kavgr=kavgr+korr;kavgi=kavgi+kori;
	  }
	kavgr=kavgr/avg;kavgi=kavgi/avg;
	fprintf(ff,"\n%9.1f,",ts);
	fprintf(ff,"%+2.5f,%+2.5f",kavgr,kavgi);
	slikar[x][y]=kavgr;
	slikai[x][y]=kavgi;
	printf("  %+1.5f %+1.5f",kavgr,kavgi);
	sound(2000);delay(200);nosound();
	}
  sound(1500);delay(200);sound(1200);delay(200);sound(1000);delay(200);nosound();
  printf("\n -------------------------------------------------");
  }

utez(size,2);
fft_2D(size);
promuckaj(size);
risisl(size,20.0);


x=abs(beri1935(0));y=abs(beri1935(1));
if ((x!=64)||(y!=64))
  {
  printf("\n Synthesizer reset during observation / recorded data may be corrupted");
  fprintf(ff,"\n Synthesizer reset during observation / recorded data may be corrupted");
  fprintf(ff,"\n status: %d %d",x,y);
  }

fclose(ff);

printf("\n Record another image? ");
gets(string);seeno = strchr(strlwr(string),yes);

  }
while (seeno);

}

//-------------------------------------------------------------

void risi_polar()                //draw correlaion on a polar graph
{
int i,x,y;
float korr,kori,kora,kork,mer;


mer=5.0;    //scaling factor
i=0;
do                                           //circular graph
  {
  if (i==0)
	{
	setviewport(50,500,10,460,1);
	clearviewport();
	setcolor(GREEN);
	for (y=45;y<226;y=y+45) circle(405,235,y);
	line(180,235,630,235);line(405,10,405,460);
	circle(80,80,70);circle(80,80,70/mer);
	setcolor(WHITE);
	gotoxy(68,2);printf("%1.3f/div %c \n",.2/mer,13);
	gotoxy(2,25);printf("\n\npress <enter> to continue");
	gotoxy(4,14);printf(" abs    ph");
	}
  vzorcenje();
  korelator(&korr,&kori);
  kora=sqrt(korr*korr+kori*kori);
  if (kori==0)
	kork=0;
  else
	kork=atan2(kori,korr)*180/3.14159;
  x=405+(int)225*korr*mer;
  y=235-(int)225*kori*mer;
  if (kora<1/mer)
	{
	line(x-5,y,x+5,y);line(x,y-5,x,y+5);
	}
  x=80+(int)70*korr;
  y=80-(int)70*kori;
  line(x-5,y,x+5,y);line(x,y-5,x,y+5);
  gotoxy(1,15+i);printf("   %1.4f %+3.1f",kora,kork);
  i=(i+1)%11;
  }
while (kbhit()==0);
getch();
}

//--------------------------------------------------------------

void pisi_korel()     //print out correlation values
{
float korr,kori,kora,kork;
int i;

clrscr();gotoxy(1,25);printf("Press <enter> to continue");
i=0;
do
  {
  gotoxy(1,i+2);
  vzorcenje();
  korelator(&korr,&kori);
  kora=sqrt(korr*korr+kori*kori);
  if (kori==0)
	kork=0;
  else
	kork=atan2(kori,korr)*180/3.14159;
  printf("Korelacija: Re= %+2.5f  Im= %+2.5f",korr,kori);
  printf("      Abs= %2.5f  Kot= %+2.5f\n",kora,kork);
  clreol();
  i=(i+1)%20;
  }
while (kbhit()==0);
getch();
}

//---------------------------------------------------------------
//Function open_datax:
//
// creates the required DataX directory structure if necessary
// opens the file with metadata for reading
// opens the DataX output file for writing
// calls pisi_meta that writes metadata to the output file

int open_datax()
{
char dir1[6],dir2[6],dir3[6],dir4[6],kjesmo[50];
struct time d_time;
struct date d_date;

char leto[6],mesec[6],dan[6],ur[6],min[6],sek[6],cas[8];

char fname[15],string[80];

//printf("\n\n ---------------------------- \n\n");

getcurdir(3,kjesmo);//printf("\n smo na %s",kjesmo);


getdate(&d_date);gettime(&d_time);
sprintf(leto, "%4d",d_date.da_year);
sprintf(mesec,"%02d",d_date.da_mon);
sprintf(dan,  "%02d",d_date.da_day);
sprintf(ur   ,"%02d",d_time.ti_hour);
sprintf(min,  "%02d",d_time.ti_min);
sprintf(sek,  "%02d",d_time.ti_sec);
strcpy(cas,ur);strcat(cas,min);strcat(cas,sek);

strcpy(dir1,leto);strcpy(dir2,mesec);strcpy(dir3,dan);strcpy(dir4,ur);

fgets(string,80,fm);      //preskoci glavo
fgets(string,80,fm);
fgets(string,4,fm);       //ident
strcpy(fname,string);strcat(fname,min);
strcat(fname,sek);strcat(fname,"I");strcat(fname,".csv");

setdisk(2);       //C:

if (chdir("\\datax"))
  {
  printf("\n Nema datax na disku %d \n",getdisk());
  perror("chdir()");
  return 1;
  }

if (chdir(dir1))
  {
  if (mkdir(dir1))
	{
	printf("\n   ne moze %s",dir1);
	return 2;
	}
  else
//	printf("\n kreiran %s",dir1);
	chdir(dir1);
  }

if (chdir(dir2))
  {
  if (mkdir(dir2))
	{
	printf("\n   ne moze %s",dir2);
	return 2;
	}
  else
//	printf("\n   kreiran %s",dir2);
	chdir(dir2);
  }

if (chdir(dir3))
  {
  if (mkdir(dir3))
	{
	printf("\n     ne moze %s",dir3);
	return 3;
	}
  else
//	printf("\n     kreiran %s",dir3);
	chdir(dir3);
  }

if (chdir(dir4))
  {
  if (mkdir(dir4))
	{
	printf("\n       ne moze %s",dir4);
	return 4;
	}
  else
//	printf("\n       kreiran %s",dir4);
	chdir(dir4);
  }

ff=fopen(fname,"w");

chdir("\\");chdir(kjesmo);
return 0;
}

//-------------------------------------------------------------
//function odrezi cuts a string at the semicolon, removes LF

void odrezi(char *string)
{
char *ppp, pp=';';

ppp = strchr(string, pp);           // poisce podpicje
string[strlen(string)-1] = '\0';    // odreze LF
if (pp) string[ppp-string] = '\0';  // odreze od podpicja dalje

}

//-------------------------------------------------------------
//function pisi_meta  writes metadata to the output DataX file

void pisi_meta(double ts, float brate, float srate, double pfrek, int rdiv1, int mdiv1, int rdiv2, int mdiv2, char sname[], char comment[], char assfile[])
{
char string[80];

fseek(fm, 0, SEEK_SET);
fgets(string,80,fm);      //preskoci glavo
fgets(string,80,fm);
fgets(string,80,fm);
fgets(string,11,fm);      //ident
fprintf(ff,"%s.SIDI11,%9.0f:DataX",string,ts);
fgets(string,80,fm);      //preskoci ostanek identa?

//global metadata follows

fprintf(ff,"\n,Source:    %s",sname);
fprintf(ff,"\n,Comment:   %s",comment);
fprintf(ff,"\n,Associated_files: %s",assfile);
fgets(string,80,fm);odrezi(string);      //HW version
fprintf(ff,"\n,HW_version: %s",string);
fprintf(ff,"\n,SW_version: 1.1.1");
fgets(string,80,fm);odrezi(string);      //data type
fprintf(ff,"\n,Data_type: %s",string);
fgets(string,80,fm);odrezi(string);      //No of channels
fprintf(ff,"\n,No_of_channels: %s",string);
fgets(string,80,fm);odrezi(string);      //No of baselines
fprintf(ff,"\n,No_of_baselines: %s",string);
fgets(string,80,fm);odrezi(string);      //reference type
fprintf(ff,"\n,Master_oscillator: %s",string);
fgets(string,80,fm);odrezi(string);      //reference error
fprintf(ff,"\n,Master_oscillator_error: %s",string);
fgets(string,80,fm);odrezi(string);      //locking
fprintf(ff,"\n,Locking: %s",string);
fgets(string,80,fm);odrezi(string);      //locking source
fprintf(ff,"\n,Locking_source: %s",string);
fgets(string,80,fm);odrezi(string);      //time stamp error
fprintf(ff,"\n,Time_stamp_error: %s",string);

//channel 1 metadata follows

fgets(string,80,fm);             //preskoci prazno vrstico
fprintf(ff,"\n,Metadata_for_channel_1");
fgets(string,80,fm);odrezi(string);      //ant type
fprintf(ff,"\n,,Ant_type: %s",string);
fgets(string,80,fm);odrezi(string);      //ant gain
fprintf(ff,"\n,,Ant_gain: %s",string);
fgets(string,80,fm);odrezi(string);      //ant polar
fprintf(ff,"\n,,Ant_polarization: %s",string);
fgets(string,80,fm);odrezi(string);      //ant mount
fprintf(ff,"\n,,Ant_mount: %s",string);
fgets(string,80,fm);odrezi(string);      //mounnt data
fprintf(ff,"\n,,Mount_data: %s",string);
fgets(string,80,fm);odrezi(string);      //ant direction
fprintf(ff,"\n,,Ant_direction: %s",string);
fgets(string,80,fm);odrezi(string);      //ant direction error
fprintf(ff,"\n,,Ant_direction error: %s",string);
fgets(string,80,fm);odrezi(string);      //position of phase center
fprintf(ff,"\n,,Position of phase center: %s",string);
fgets(string,80,fm);odrezi(string);      //position error
fprintf(ff,"\n,,Position error: %s",string);
fgets(string,80,fm);odrezi(string);      //electric delay
fprintf(ff,"\n,,Electric_delay: %s",string);
fgets(string,80,fm);odrezi(string);      //electric delay error
fprintf(ff,"\n,,Electric_delay_error: %s",string);
fgets(string,80,fm);odrezi(string);      //phase correction
fprintf(ff,"\n,,Phase_correction: %s",string);
fgets(string,80,fm);odrezi(string);      //phase error
fprintf(ff,"\n,,Phase error: %s",string);
fgets(string,80,fm);odrezi(string);      //system noise temperature
fprintf(ff,"\n,,System_noise_temperature: %s",string);
fgets(string,80,fm);odrezi(string);      //noise temperature error
fprintf(ff,"\n,,Noise_temperature_error: %s",string);
fgets(string,80,fm);odrezi(string);      //channel type
fprintf(ff,"\n,,Channel_type: %s",string);
fgets(string,80,fm);odrezi(string);      //sampling
fprintf(ff,"\n,,Sampling: %s",string);
fgets(string,80,fm);odrezi(string);      //A/D no of bits
fprintf(ff,"\n,,A/D_no_of_bits: %s",string);
fgets(string,80,fm);odrezi(string);      //analog filter
fprintf(ff,"\n,,Analog_filter: %s",string);
fgets(string,80,fm);odrezi(string);      //data packing method
fprintf(ff,"\n,,Data_packing: %s",string);
fprintf(ff,"\n,,Sample_rate: %8.0f",srate);
fprintf(ff,"\n,,Sample_rate_ratios: 0,0");
fprintf(ff,"\n,,Data_rate: %2.3f",brate);
fprintf(ff,"\n,,Samples_per_datapoint: 131072");
fprintf(ff,"\n,,Frequency: %4.6f",pfrek);
fprintf(ff,"\n,,Frequency_ratios: %d,%d",rdiv1,mdiv1);


//channel 2 metadata follows

fgets(string,80,fm);             //preskoci prazno vrstico
fprintf(ff,"\n,Metadata_for_channel_2");
fgets(string,80,fm);odrezi(string);      //ant type
fprintf(ff,"\n,,Ant_type: %s",string);
fgets(string,80,fm);odrezi(string);      //ant gain
fprintf(ff,"\n,,Ant_gain: %s",string);
fgets(string,80,fm);odrezi(string);      //ant polar
fprintf(ff,"\n,,Ant_polarization: %s",string);
fgets(string,80,fm);odrezi(string);      //ant mount
fprintf(ff,"\n,,Ant_mount: %s",string);
fgets(string,80,fm);odrezi(string);      //mounnt data
fprintf(ff,"\n,,Mount_data: %s",string);
fgets(string,80,fm);odrezi(string);      //ant direction
fprintf(ff,"\n,,Ant_direction: %s",string);
fgets(string,80,fm);odrezi(string);      //ant direction error
fprintf(ff,"\n,,Ant_direction error: %s",string);
fgets(string,80,fm);odrezi(string);      //position of phase center
fprintf(ff,"\n,,Position of phase center: %s",string);
fgets(string,80,fm);odrezi(string);      //position error
fprintf(ff,"\n,,Position error: %s",string);
fgets(string,80,fm);odrezi(string);      //electric delay
fprintf(ff,"\n,,Electric_delay: %s",string);
fgets(string,80,fm);odrezi(string);      //electric delay error
fprintf(ff,"\n,,Electric_delay_error: %s",string);
fgets(string,80,fm);odrezi(string);      //phase correction
fprintf(ff,"\n,,Phase_correction: %s",string);
fgets(string,80,fm);odrezi(string);      //phase error
fprintf(ff,"\n,,Phase error: %s",string);
fgets(string,80,fm);odrezi(string);      //system noise temperature
fprintf(ff,"\n,,System_noise_temperature: %s",string);
fgets(string,80,fm);odrezi(string);      //noise temperature error
fprintf(ff,"\n,,Noise_temperature_error: %s",string);
fgets(string,80,fm);odrezi(string);      //channel type
fprintf(ff,"\n,,Channel_type: %s",string);
fgets(string,80,fm);odrezi(string);      //channel type
fprintf(ff,"\n,,Recorded_signal: %s",string);
fgets(string,80,fm);odrezi(string);      //sampling
fprintf(ff,"\n,,Sampling: %s",string);
fgets(string,80,fm);odrezi(string);      //A/D no of bits
fprintf(ff,"\n,,A/D_no_of_bits: %s",string);
fgets(string,80,fm);odrezi(string);      //analog filter
fprintf(ff,"\n,,Analog_filter: %s",string);
fgets(string,80,fm);odrezi(string);      //data packing method
fprintf(ff,"\n,,Data_packing: %s",string);
fprintf(ff,"\n,,Sample_rate: %8.0f",srate);
fprintf(ff,"\n,,Sample_rate_ratios: 0,0");
fprintf(ff,"\n,,Data_rate: %2.3f",brate);
fprintf(ff,"\n,,Samples_per_datapoint: 131072");
fprintf(ff,"\n,,Frequency: %4.6f",pfrek);
fprintf(ff,"\n,,Frequency_ratios: %d,%d",rdiv2,mdiv2);

fprintf(ff,"\n,Correlation_values_for_baseline_1");
fprintf(ff,"\nTime,Real,Imag");

fclose(fm);

}

//-----------------------------------------------------------------
//asks user for frequency and sets the synthesizers

void set_freq(double *frek, double *pfrek, int *rdiv1, int *mdiv1, int *rdiv2, int *mdiv2)
{
double frk,pfrk;
int r1,m1,r2,m2;
int st1,st2;

printf("\n\n Enter frequency! [MHz] ");scanf("%lf",&frk);fflush(stdin);
i2creset();
st1=sl1935(0,frk,&pfrk,&r1,&m1);st2=sl1935(1,frk,&pfrk,&r2,&m2);
printf("\n");
if ((st1&15)!=0)
  {
  printf("\nChannel 1 synthesizer not responding!");
  }
if ((st2&15)!=0)
  {
  printf("\nChannel 2 synthesizer not responding!");
  }

if ((st1&64)==0)
  {
  printf("\nChannel 1 synthesizer didn't lock at %lf MHz",frk);
  }
if ((st2&64)==0)
  {
  printf("\nChannel 2 synthesizer didn't lock at %lf MHz",frk);
  }

*frek=frk;*pfrek=pfrk;
*rdiv1=r1;*mdiv1=m1;
*rdiv2=r2;*mdiv2=m2;
}

//---------------------------------------------------------------
//perform some basic diagnostics

void tests(float *brate, float *srate)
{
long j;
float br,sr;

checkspeed(&br,&sr);
clrscr();

gotoxy(10,2);printf("Sample rate: %9.0f/sec     Block rate: %2.2f/sec",sr,br);
gotoxy(1,4);printf("           DC offsets                Quadrature errors  Input samples");
gotoxy(1,5);printf("   1I       1Q       2I       2Q         ch1     ch2         HEX");
gotoxy(10,24);printf("Press <enter> to continue");

j=0;
do
  {
  gotoxy(1,j+6);
  vzorcenje();checkdc();checkquad();printf("      %6X",sig[555]);printf("\n");
  clreol();
  j=(j+1)%16;
  }
while (kbhit()==0);
getch();

*brate=br;*srate=sr;
}

//--------------------------------------------------------------

void fft(unsigned int n, float xr[], float xi[], float wr[], float wi[], int bi[])
{
unsigned int b,dp,a,t,p;
float xwr,xwi;

for (b=0;b<n;b++)
  if (bi[b]>b)
	{xwr=xr[b];xwi=xi[b];
	xr[b]=xr[bi[b]];xi[b]=xi[bi[b]];
	xr[bi[b]]=xwr;xi[bi[b]]=xwi;}

b=1;dp=n>>1;
do
  {a=0;
  do
	{t=a+b;p=0;
	do
	  {xwr=xr[a+b]*wr[p]-xi[a+b]*wi[p];
	  xwi=xr[a+b]*wi[p]+xi[a+b]*wr[p];
	  xr[a+b]=xr[a]-xwr;xi[a+b]=xi[a]-xwi;
	  xr[a]=xr[a]+xwr;xi[a]=xi[a]+xwi;
	  p=p+dp;a=a++;}
	while (a<t);
	a=a+b;}
  while (a<n);
  dp=dp>>1;b=b<<1;}
while (b<n);
}


//--------------------------------------------------------------------

void fftmojprip(unsigned int n,float wr[], float wi[], int bi[])
{
unsigned int i,j,k,l,m;
m=(int)(log((double)n)/log(1.9999999999999));
for (l=0;l<n;l++)    //bit reverse adrese
  {i=0;  j=l<<1;
  for (k=0;k<m;k++)
	i=(j=j>>1)%2+(i=i<<1);
  bi[l]=i;}

for (l=0;l<n;l++)    //sin/cos tabela
  {wr[l]=cos((float)l/(float)n*2*PI);
  wi[l]=-sin((float)l/(float)n*2*PI);}
}

//-----------------------------------------

void setpal()
{
/*
setrgbpalette(1,0,0,16);
setrgbpalette(2,0,0,32);
setrgbpalette(3,0,0,48);
setrgbpalette(4,0,0,63);
setrgbpalette(5,0,16,63);
setrgbpalette(20,0,32,32);
setrgbpalette(7,0,48,16);
setrgbpalette(56,0,63,0);
setrgbpalette(57,48,63,0);
setrgbpalette(58,63,63,0);
setrgbpalette(59,63,32,0);
setrgbpalette(60,63,0,0);
setrgbpalette(61,63,16,32);
setrgbpalette(62,63,32,48);
setrgbpalette(63,63,63,63);
*/
setrgbpalette(1,4,4,4);
setrgbpalette(2,8,8,8);
setrgbpalette(3,12,12,12);
setrgbpalette(4,16,16,16);
setrgbpalette(5,20,20,20);
setrgbpalette(20,24,24,24);
setrgbpalette(7,28,28,28);
setrgbpalette(56,32,32,32);
setrgbpalette(57,36,36,36);
setrgbpalette(58,40,40,40);
setrgbpalette(59,44,44,44);
setrgbpalette(60,48,48,48);
setrgbpalette(61,52,52,52);
setrgbpalette(62,56,56,56);
setrgbpalette(63,63,63,63);

}

//-------------------------------------------------------------

void risisl(int vm, float din)     //din = dinamika dB
{
int gdriver = DETECT, gmode;
float max,min,sa,sl;
int i,j;

initgraph(&gdriver, &gmode, "c:\\bc\\bgi");

max=0;min=1e12;
for (i=0;i<vm;i++)
  {
  for (j=0;j<vm;j++)
	{
	sa=sqrt(slikar[i][j]*slikar[i][j]+slikai[i][j]*slikai[i][j]+1e-12);
	if (sa>max) max=sa;
	if (sa<min) min=sa;
	}
  }

setcolor(WHITE);moveto(40,40);lineto(60+20*vm,40);lineto(60+20*vm,60+20*vm);
lineto(40,60+20*vm);lineto(40,40);

setpal();

for (i=0;i<vm;i++)
  {
  for (j=0;j<vm;j++)
	{
	sa=sqrt(slikar[i][j]*slikar[i][j]+slikai[i][j]*slikai[i][j]+1e-12);
//	sl=15.0*(sa-min)/(max-min);
	sl=15.0+10.0*log10(sa/max)*15.0/din;
	if (sl<0.0) sl=0.0;
	setfillstyle(SOLID_FILL,(int)sl);
	bar(50+20*i,50+20*j,70+20*i,70+20*j);
	}
  }

for (i=0;i<16;i++)
  {
  setfillstyle(SOLID_FILL,i);
  bar(50+20*i,450,70+20*i,470);
  }

getch();
closegraph();

}

//----------------------------------------------------------------

void fft_2D(int vm)
{
int i,j,n;
n=vm;

fftmojprip(n,wr,wi,bi);

for (j=0;j<vm;j++)         //FFT po vrsticah
  {
  for (i=0;i<vm;i++)
	{
	xr[i]=slikar[i][j];
	xi[i]=slikai[i][j];
	}
  fft(n,xr,xi,wr,wi,bi);
  for (i=0;i<vm;i++)
	{
	slikar[i][j]=xr[i];
	slikai[i][j]=xi[i];
	}
  }

for (i=0;i<vm;i++)         //FFT po stolpcih
  {
  for (j=0;j<vm;j++)
	{
	xr[j]=slikar[i][j];
	xi[j]=slikai[i][j];
	}
  fft(n,xr,xi,wr,wi,bi);
  for (j=0;j<vm;j++)
	{
	slikar[i][j]=xr[j];
	slikai[i][j]=xi[j];
	}
  }

}

//-------------------------------------------------------

void utez(int vm, int tip)
{
int i,j;
float u;
float x,y;

for (i=0;i<vm;i++)
  {
  for (j=0;j<vm;j++)
	{
	x=((float)i+0.5)/(float)vm;
	y=((float)j+0.5)/(float)vm;
	switch (tip)
	  {
	  case 0:
		u=1;break;
	  case 1:
		u=sin(x*PI)*sin(y*PI);break;                             //polperioda
	  case 2:
		u=(cos((x-0.5)*2*PI)+1)*(cos((y-0.5)*2*PI)+1)/4;break;   //raised cos
	  default:
		u=0;break;
	  }
	slikar[i][j]=slikar[i][j]*u;
	slikai[i][j]=slikai[i][j]*u;
	}
  }

}

//-------------------------------------------------------

void promuckaj(int vm)
{
int i,j,vmp;
float zr,zi;

vmp=vm/2;
for (j=0;j<vm;j++)           //po vrsticah
  {
  for (i=0;i<vmp;i++)
	{
	zr=slikar[j][i];zi=slikai[j][i];
	slikar[j][i]=slikar[j][i+vmp];slikai[j][i]=slikai[j][i+vmp];
	slikar[j][i+vmp]=zr;slikai[j][i+vmp]=zi;
	}
  }

for (j=0;j<vm;j++)           //po stolpcih
  {
  for (i=0;i<vmp;i++)
	{
	zr=slikar[i][j];zi=slikai[i][j];
	slikar[i][j]=slikar[i+vmp][j];slikai[i][j]=slikai[i+vmp][j];
	slikar[i+vmp][j]=zr;slikai[i+vmp][j]=zi;
	}
  }

}

//-----------------------------------------------------------

void cakitipko()        //pocaka na stik serial port 9 pin pina 7 in 8
{                       //ali tipko na tastaturi

int mcon,msta;
char mst;

mcon=0x3FC;msta=0x3FE;     //3FC,3FE   ali   2FC,2FE

outportb(mcon,15);         //da so izhodi +

if (kbhit()) getch();             //ce je kaksen ostal od prej

while ( ((inportb(msta)&16)!=16)&&(kbhit()==0) )
  {
  }

}

//********************************************************************

int main()
{

float brate,srate;
char string[80],yes='y',*imaname;
double frek,pfrek;
int gdriver = DETECT, gmode;
int rdiv1,rdiv2,mdiv1,mdiv2;



mkcortab(0,ct0real,ct0imag);    //make lookup tables for the correlator
mkdctab(dctab);

clrscr();printf("\n\n         ----- SIDI 1.1 IMAGER starting ------ \n");

set_freq(&frek, &pfrek, &rdiv1, &mdiv1, &rdiv2, &mdiv2);

tests(&brate,&srate);

irqmask=inportb(0x21);outportb(0x21,irqmask+1);      //stop timer interrupt

//pisi_korel();              //print out some correlation values

initgraph(&gdriver, &gmode, "c:\\bc\\bgi");

risi_polar();		//cicular graph

closegraph();

clrscr();gotoxy(1,1);

printf("\n Record data to file? ");
gets(string);imaname = strchr(strlwr(string),yes);

if (imaname)
  {
  risi_franze(brate,srate,pfrek,rdiv1,mdiv1,rdiv2,mdiv2);      //draw fringes
  }

outportb(0x21,irqmask);      //resume timer interrupt
nastavi_uro();               //restore time/date from HW clock


return 0;
}