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

#include "tfitshdr.h"
#include "inifile.h"
#define CHANNELS 120

char *btype[]={"date","time","object",
"freq","cdelt","crpix","tqsun","theta","kflux",
"solar_b","solar_p","solar_r","azimuth","altitude",
"fiv","fcalibr"};

char *bform[]={"16A","16A","16A","E","E","E","E","E",
"E","E","E","E","E","E","I","I"};

typedef struct
{
char date[16],time[16],object[16];
float freq,cdelt,crpix,tqsun,theta,kflux;
float solar_b,solar_p,solar_r,azimuth,altitude;
short fiv,fcalibr;
} SbinFits;

typedef struct
{
 short size,regOn;
 unsigned int startDT, timetik;
 unsigned char modAuto, modManual, ngOn, ngLeft;
 unsigned char ngRight, attMask1, attMask2, attMask3;
 int aver_max_value;
 int data[240];
} TArmData;

float freq[]={
758.8, 801.8, 846, 891.5, 938.3, 986.4, 1035.9, 1086.9, 1139.4, 1193.4, 1249, 1306.2, 1365, 1425.5,
1487.8, 1551.9, 1617.9, 1685.9, 1755.9, 1827.9, 1901.9, 1978, 2056.3, 2136.8, 2219.6, 2304.8, 2392.5,
2482.8, 2570, 2670, 2770, 2890, 3000, 3180, 3370, 3570, 3750, 4000, 4200, 4400, 4500, 4700, 4900, 5100,
5300, 5500, 5700, 6000, 6150, 6300, 6650, 6800, 6950, 7100, 7300, 7600, 7750, 7900, 8100, 8250, 8400,
8700, 8900, 9050, 9200, 9350, 9500, 9700, 10000, 10150, 10300, 10500, 10600, 10800, 10900, 11300,
11400, 11600, 11700, 11900, 12200, 12400, 12500, 12700, 12900, 13050, 13350, 13500, 13700, 14000,
14150, 14300, 14500, 14600, 15000, 15300, 15450, 15600, 15750, 15900, 16100, 16250, 16400, 16500,
16900, 17050, 17200, 17350, 17500, 17700, 17800, 18200, 0, 0, 0, 0, 0, 0, 0, 0 
};

char* str2upper(char *x){
 static char name[64];
 int i,j=strlen(x);
 if(j>63) j=63;
 for(i=0;i<j;i++)name[i]=toupper(x[i]);name[j]=0;
 return name;
}

//Tqsun = -0,0042x3 + 0,1991x2 - 6,9102x + 6117,5 where x=freq in GGz
float Tqsun(float x)
{
 return -0.0042*x*x*x + 0.1991*x*x - 6.9102*x + 6117.5;
}

//Flux_koef = 0,0044x + 0,0308  where x=freq in GGz
float fKoef(float x)
{
 return 0.0044*x + 0.0308;
}

//theta=float(0.2+0.94*3.0e2/freq[GGz]);
float Theta(float x)
{
if(x) return 0.2+0.94*300/x; else return 1;
}

void
ScanMove(int *dat, float move, int nn)
{
 int j=0,i,x;
 float z;
    if(move==0) return;
    x=(int)move;
    if((move-(float)x)==0) j=1;
    if(x<0)
      for(i=nn+x-1;i>=0;i--) dat[i-x]=dat[i];
    else
      if(x!=0) for(i=x;i<nn;i++) dat[i-x]=dat[i];
 if(j) return;
  z=fabs(move-(float)x);
    if(move<0)
      for(i=nn-2;i>=0;i--) dat[i+1]+=(dat[i]-dat[i+1])*z;
    else
      for(i=1;i<nn;i++) dat[i-1]+=(dat[i]-dat[i-1])*z;
}

int main(int argc, char* argv[])
{
TArmData *armData;
TIniFile ini;
int i, j, k, l, m, d, ArmDataNum, ArmChannelNum, smooth=0, armSmooth, fdebug=0, eoftest=0;
float lowshift, hishift, shiftfrq, fshift;
int lostBlock[3][1000], lostBlockNum, lostpix, naxis1=0, *rawData, *fitsData, snum;
unsigned int startobs, stopobs, size, hdrsize, ii;
unsigned char *data;
char *c, stmp[256], date[32], time[32], obj[64], FreqFile[256], iniDirFile[256], cc;
char dataFileMask[256], dataDir[256], fitsDir[256], sunCoef[256];
float z, armDT, azimuth=0, altitude=0, sol_dec=0, sol_ra=0, solar_r=900, solar_p=0, solar_b=0, sol_valh=0;
float frqFlag[2][CHANNELS], crpix, tbsun, kflux, theta;
double cdelt1, sum;
TFitsHdr fits,fraw;

FreqFile[0]=0; iniDirFile[0]=0; sunCoef[0]=0;
QDateTime DateTime;
DateTime=QDateTime::currentDateTime();
{
 QString s("yyyy-MM-dd_hh:mm:ss");
 s=DateTime.toString(s);
 QByteArray ss;
 ss.append(s);
 c=(char *)ss.constData();
 printf("arm2fits starts at %s\n",c);
}
if(argc<2){printf("Error - no arguments\n"); return -1;}
// if(argc<2)sprintf(dataFileMask,"moon0"); else
sprintf(dataFileMask,"%s",argv[1]);
sprintf(stmp,"%s.ini",argv[0]);
printf(" iniFile=%s\n dataFileMask=%s\n",stmp,dataFileMask);
if(ini.LoadIniFile(stmp))
{
 fdebug=ini.GetParam("debug",(int)0);
 smooth=ini.GetParam("smooth",(int)1);
 eoftest=ini.GetParam("eoftest",(int)1);
 lowshift=ini.GetParam("lowshift",(float)0);
 hishift=ini.GetParam("hishift",(float)0);
 shiftfrq=ini.GetParam("shiftfrq",(float)3.1);
 if(smooth<1) smooth=1;
 ini.GetParam("FreqFile",FreqFile);
 ini.GetParam("sunCoef",sunCoef);
 ini.GetParam("iniDirFile",iniDirFile);
 if(iniDirFile[0]<33)strcpy(iniDirFile,stmp);
 if(fdebug) printf(" FreqFile=%s\n iniDirFile=%s\n",FreqFile,iniDirFile);
 ini.Clear();
} else {printf("Error - can not open program INI file\n"); return -1;} 
if(ini.LoadIniFile(stmp))
{
 dataDir[0]=0; ini.GetParam("dataDir",dataDir);
 fitsDir[0]=0; ini.GetParam("fitsDir",fitsDir);
 if(fdebug) printf("dataDir=%s\n fitsDir=%s\n",dataDir,fitsDir);
 ini.Clear();
} else {printf("Error - can not open data INI file\n"); return -1;} 

if(FreqFile[0]<33) 
{
 for(i=0;i<CHANNELS;i++){frqFlag[1][i]=freq[i];if(freq[i])frqFlag[0][i]=1; else frqFlag[0][i]=0;}
 printf("Frequensy configuration is taken from internal memory\n");
}
else
{
 ini.Clear();
 float l[8];
 j=ini.LoadIniFile(FreqFile);
 if(j==0){printf("Error - can not open channels configuration file %s\n",FreqFile); return -1;}
 for(i=0;i<CHANNELS;i++)
 {
  z=i+1;
  j=ini.GetParam(z,l);
  if(j==3){frqFlag[1][i]=l[0];frqFlag[0][i]=l[1];}
  else {frqFlag[1][i]=0;frqFlag[0][i]=0;}
  if(fdebug) printf("freq=%f flag=%f\n",frqFlag[1][i],frqFlag[0][i]);
 }
 printf("Frequensy configuration is taken from %s\n",FreqFile);
 ini.Clear();
}
for(i=0,j=0;i<CHANNELS;i++)if(frqFlag[0][i])j++;
ArmChannelNum=j;
printf("Frequensy configuration has %d channels\n",ArmChannelNum);
//-- opening RAW data file ---
sprintf(stmp,"%s%s.raw",dataDir,dataFileMask);
QFile raw(stmp);
if(!raw.open(QIODevice::ReadOnly)){printf("Error - can not open <%s>\n", stmp); return -1;}
else printf("RAW data file <%s> is opened for reading\n", stmp);
size=raw.size(); 
data=raw.map(0L,size);
if(data==0) {printf("Error - can not map raw file to memory\n"); return -1;}
hdrsize=fraw.LoadFromMem((char*)data);
if(hdrsize<2880){printf("Error - can not read mapped memory\n"); return -1;}
stmp[0]=0;
fraw.GetParam("extname",stmp);
if(strstr(stmp,"Row_reg")!=stmp){printf("Error - raw file format is wrong\n"); return -1;}
ArmDataNum=(size-hdrsize)/sizeof(TArmData);
printf("RAW data file has %d data blocks\n", (int)ArmDataNum);
if(eoftest) if(size-hdrsize-ArmDataNum*sizeof(TArmData)){printf("Error - raw file size is wrong\n"); return -1;}
armData=(TArmData *)(data+hdrsize);
startobs=armData[0].startDT;
armSmooth=armData[0].aver_max_value;
if(armSmooth<1) armSmooth=1;
stopobs=armData[ArmDataNum-1].startDT;
//-- opening output FITS file ---
sprintf(stmp,"%s%s.fits",fitsDir,dataFileMask);
QFile out(stmp);
if(!out.open(QIODevice::WriteOnly)){printf("Error - can not open <%s>\n", stmp); return -1;}
else printf("FITS file <%s> is opened for writing\n", stmp);

//--- lost data block check ----
lostBlock[0][0]=0;
for(j=0,i=1;i<(ArmDataNum-1);i++)
{
 k=armData[i].timetik-armData[i-1].timetik;
 if(k>1)
 {
  lostBlock[1][j]=i;lostBlock[2][j]=k-1;  
  if(fdebug) printf("data frame %d was lost\n",i);
  j++;
  lostBlock[0][j]=i;
 }
 if(j>=999) {printf("Error - More than 1000 data blocks have been lost\n"); return -1;} 
}
lostBlock[1][j]=ArmDataNum; lostBlock[2][j]=0;
lostBlockNum=j;
lostpix=0;
if(lostBlockNum)for(i=0;i<lostBlockNum;i++)lostpix+=lostBlock[2][i];
naxis1=ArmDataNum+lostpix;
if(lostBlockNum) printf("%d data blocks (%d pixels) have been lost\n",lostBlockNum,lostpix);
else printf("No lost blocks in data file\n");
if(fdebug)
{
 if(lostBlockNum) 
 for(i=0;i<=lostBlockNum;i++) printf("%d range=%d-%d delta=%d\n",i,lostBlock[0][i],lostBlock[1][i],lostBlock[2][i]);
 printf("ArmChannelNum=%d ArmDataNum=%d\n",ArmChannelNum,ArmDataNum);
 printf("tik0=%d tik1=%d tik2=%d\n",armData[0].timetik,armData[ArmDataNum-2].timetik,armData[ArmDataNum-1].timetik);
 printf("startobs=%d stopobs=%d\n",startobs,stopobs);
 i=stopobs-startobs;z=i;
 printf("delta=%d cdelt=%f\n",i, z/ArmDataNum);
}

//--- FITS header creation ----
printf("\nCreating output FITS file Header...\n");
fraw.GetParam("date-obs",date);
fraw.GetParam("time-obs",time);
if(fraw.GetParam("armdt",armDT)) armDT=0.007;
fraw.GetParam("object",obj);
fraw.GetParam("azimuth",azimuth);   
fraw.GetParam("altitude",altitude);
fraw.GetParam("sol_dec",sol_dec);
fraw.GetParam("sol_ra",sol_ra);
fraw.GetParam("solar_r",solar_r);
fraw.GetParam("solar_p",solar_p);
fraw.GetParam("solar_b",solar_b);
fraw.GetParam("sol_valh",sol_valh);
//--- cdelt calculation --- 
cdelt1=tan(sol_dec*0.017453293)*tan(azimuth*0.017453293);
cdelt1=sqrt(1-cdelt1*cdelt1);
cdelt1*=15.0*cos(sol_dec*0.017453293)*(1.0-sol_valh/3600.0);
k=stopobs-startobs;
if(k>0) cdelt1*=k/(double)naxis1; else cdelt1*=armDT*armSmooth;
if(smooth>1) cdelt1*=smooth;

//--- crpix calculation ---
if(smooth>1) i=naxis1/smooth; else i=naxis1;
sprintf(stmp,"%s_%s",date,time);
{
 QString s(stmp);
 QString s1("yyyy/MM/dd_HH:mm:ss.zzz");
 DateTime=QDateTime::fromString(s,s1); 
 ii=DateTime.toTime_t(); m=ii%3600; j=startobs%3600; if(m<j)m+=3600;
 sum=m-j;
 if(fdebug)printf("%s m=%d sum=%d\n",stmp,m,(int)sum);
}
if(k>0) z=double(k)/i; else {z=armDT*armSmooth; if(smooth>1) z*=smooth;}
sum/=z;
crpix=i-sum-1;
if(fdebug)printf("cdelt=%f deltax=%f crpix=%f\n",z,(float)sum,crpix);
if (crpix <0 ||crpix > i) crpix=i/2.0;  
//--- FITS Header creation ----
fits.makeHdr("");
fits.PutParam("bitpix",32);
fits.PutParam("naxis",3);
fits.PutParam("naxis1",i);
fits.PutParam("naxis2",2);
fits.PutParam("naxis3",ArmChannelNum);
fits.PutParam("dataname","newpas");
fits.PutParam("TELESCOP","RATAN-600");
fits.PutParam("origin","PAS-120");
fits.PutParam("date-obs",date);
fits.PutParam("time-obs",time);
fits.PutParam("startobs",(int)startobs);
fits.PutParam("stopobs",(int)stopobs);
fits.PutParam("nsamples",ArmDataNum,"NAXIS1 = (NSAMPLES + lostpixs)/SMOOTH");
fits.PutParam("lostpixs",lostpix);
fits.PutParam("smooth", smooth);
fits.PutParam("armsmoot",armSmooth);
fits.PutParam("armdt",armDT);
fits.PutParam("crpix1",crpix);
fits.PutParam("cdelt1",cdelt1,"arcsec");
fits.PutParam("flag_iv",1,"1- data is RL; 0- data is IV");
fits.PutParam("Calibr",0);
fits.PutComment("*** Object parameters ***");
fits.PutParam("object",obj);
fits.PutParam("azimuth",azimuth);   
fits.PutParam("altitude",altitude);
fits.PutParam("sol_dec",sol_dec);
fits.PutParam("sol_ra",sol_ra);
fits.PutParam("solar_r",solar_r,"arcsec");
fits.PutParam("solar_p",solar_p);
fits.PutParam("solar_b",solar_b);
fits.PutParam("sol_valh",sol_valh);
fits.PutComment("*** Frequencies ***");
for(i=0,j=0;i<CHANNELS;i++)
 if(frqFlag[0][i])
 {
   j++;
   sprintf(stmp,"FREQ%03d",j);
   fits.PutParam(stmp,float(frqFlag[1][i]/1000),"GHz");
 } 
c=fits.WriteToMem(i);
if(c)
{
out.write(c,i);
delete[] c;
} 
else {printf("Error - can not write FITS header to file\n"); return -1;}
printf("FITS data dimention is %d x 2 x %d (smooth=%d)\n", int(naxis1/smooth),ArmChannelNum,smooth);
//--- Writing arm data to FITS ---
printf("Writing arm data to FITS...\n");
rawData= new int[naxis1]; 
if(smooth >1) {snum=(int)(naxis1/smooth); fitsData=new int[snum];} else {fitsData=rawData; snum=naxis1;}
for(i=0;i<CHANNELS;i++)
 if(frqFlag[0][i]) for(j=0;j<2;j++)
 {
   if(lostBlockNum)
   {
    d=0;
    for(k=0;k<=lostBlockNum;k++)
	 {
          for(l=lostBlock[0][k];l<lostBlock[1][k];l++,d++) rawData[d]=armData[l].data[2*i+j];
	  m=rawData[d-1];
          for(l=0;l<lostBlock[2][k];l++,d++) rawData[d]=m;
	 }
	if(d<naxis1) printf("%d < %d(naxis1)\n",d,naxis1); 
   }
   else for(k=0;k<ArmDataNum;k++) rawData[k]=armData[k].data[2*i+j];
   // ---here should be shifting procedure----
   if(frqFlag[1][i]>shiftfrq) fshift=hishift; else fshift=lowshift;
   if(j) fshift*=-1; fshift/=armSmooth;
   ScanMove(rawData,fshift,naxis1);
   if(fdebug) printf("freq=%f shift=%f\n",frqFlag[1][i],fshift);
   //----------------------------------
   if(smooth>1)
    {
        for(k=0;k<snum;k++){for(d=0,sum=0;d<smooth;d++) sum+=rawData[k*smooth+d]; sum/=smooth; fitsData[k]=(int)sum;}
	}
   for(k=0;k<(snum/2);k++){d=fitsData[k];fitsData[k]=fitsData[snum-1-k];fitsData[snum-1-k]=d;}
   c=(char*)fitsData; 
   for(k=0;k<snum;k++)
    {
    cc=c[k*4];c[k*4]=c[k*4+3];c[k*4+3]=cc;
    cc=c[k*4+1];c[k*4+1]=c[k*4+2];c[k*4+2]=cc;
	}
   out.write(c,snum*4);
 }
d=ArmChannelNum*2*snum*4;
if(d%2880){k=2880-d%2880; c=new char[k]; memset(c,0,k); out.write(c,k);delete[] c;} 
if(fdebug)printf("d=%d k=%d\n",d%2880, k);
if(smooth >1) delete[] fitsData; 
delete[] rawData; 

//--- create Binary table extension --- 
fits.Clear();
fits.makeHdr("BINTABLE");
printf("\nCreating BINTABLE extension...\n");
fits.PutParam("EXTNAME","Scan_params");
fits.PutParam("naxis",2);
i=sizeof(SbinFits);
fits.PutParam("naxis1",i,"Number of bytes per row");
fits.PutParam("naxis2",ArmChannelNum,"Number of rows");
fits.PutParam("pcount",0,"Normally 0 (no varying arrays)");
fits.PutParam("gcount",1);
fits.PutParam("tfields",16,"Number of columns in table");
fits.PutComment("*** Column names ***");
for(i=0;i<16;i++){
 sprintf(stmp,"TTYPE%d",i+1);
 fits.PutParam(stmp,btype[i]);
}
fits.PutComment("*** Column formats ***");
for(i=0;i<16;i++){
 sprintf(stmp,"TFORM%d",i+1);
 fits.PutParam(stmp,bform[i]);
}
c=fits.WriteToMem(i);
if(c)
{
out.write(c,i);
delete[] c;
} 
else {printf("Error - can not write BINFITS header to file\n"); return -1;}

//--- writing Binary table extension --- 
SbinFits bstr;
c=(char*)&bstr;
d=sizeof(SbinFits); memset(&bstr,0,d);
i=strlen(date);if(i>15) i=15; memmove(bstr.date,date,i);
i=strlen(time);if(i>15) i=15; memmove(bstr.time,time,i);
i=strlen(obj);if(i>15) i=15; memmove(bstr.object,obj,i);
bstr.solar_p=fitsnum(float(solar_p));
bstr.solar_b=fitsnum(float(solar_b));
bstr.solar_r=fitsnum(float(solar_r));
bstr.cdelt=fitsnum(float(cdelt1));
bstr.crpix=fitsnum(crpix);
bstr.azimuth=fitsnum(float(azimuth));
bstr.altitude=fitsnum(altitude);
bstr.fiv=fitsnum(short(1));
if(strstr(str2upper(obj),"MOON"))tbsun=220.0; else tbsun=6000;
if(sunCoef[0])
 { 
 j=ini.LoadIniFile(sunCoef);
 if(j==0) printf("Error - can not open %s\n Tqsun, FluxKoef, Theta will be calculated analytically\n",sunCoef);
 else printf("Tqsun, FluxKoef, Theta are taken from %s\n", sunCoef);
 } 
else {j=0; printf("Tqsun, FluxKoef, Theta are calculated analytically\n");}
float sunC[8]; for(i=0;i<8;i++) sunC[i]=0;
for(i=0;i<CHANNELS;i++) if(frqFlag[0][i])
{
 z=frqFlag[1][i]/1000;
 if(j)
  {
  float ftmp;
  int l0=ini.GetLineInterpolation(z,sunC); 
  ftmp=0; if(l0>0) ftmp=sunC[0]; if(ftmp==0) ftmp=Tqsun(z); if(tbsun>300) tbsun=ftmp;
  ftmp=0; if(l0>1) ftmp=sunC[1]; if(ftmp==0) ftmp=fKoef(z); kflux=ftmp;
  ftmp=0; if(l0>2) ftmp=sunC[2]; if(ftmp==0) ftmp=Theta(z); theta=ftmp;
  }
 else
  {
  if(tbsun>300) tbsun=Tqsun(z);
  kflux=fKoef(z); 
  theta=Theta(z);
  }
 bstr.freq=fitsnum(z);
 bstr.kflux=fitsnum(kflux);
 bstr.theta=fitsnum(theta);
 bstr.tqsun=fitsnum(tbsun);
 if(fdebug) printf("\tfreq=%f fluxkoef=%f Tqsun=%f theta=%f\n",z,kflux,tbsun,theta);
 out.write(c,d);
}
d*=ArmChannelNum;
if(d%2880){k=2880-d%2880; c=new char[k]; memset(c,0,k); out.write(c,k);delete[] c;} 
if(fdebug)printf("d=%d k=%d\n",d%2880, k);
printf("FITS file is created\n");
return 0;
}