#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include "kg.h"
#include "SpiceUsr.h"
#include "SpiceZfc.h"

#define X0 (64+1200+64)
#define X1 (64)
#define X2 (64+1200)
#define Y0 (24+600+20+24)
#define Y1 (24)
#define Y2 (24+600)
#define Y3 (24+600+20)

#define GLAT ( 37.700)
#define GLON (140.667)

#define NS (900)
#define P1 (-110)
#define P2 ( 110)

#define NL (-300)
//#define SS (100)

static double WF2[4096];

void fft4096_init(void) {
  int i;
  for (i=0;i<4096;i++) WF2[i]=0.5-0.5*cos((2*M_PI/4096)*i);
}

int rv4096(int x) {
  int i,y;
  y=0;
  for (i=0;i<12;i++) {
    y*=2;
    if ((x>>i)&1) y+=1;
  }
  return y;
}

void fft4096(double *sw,double *PWR) {
  int i,j,k,a;
  double xr0,xi0,xr1,xi1,x0,x1,x2,x3,cth,sth;
  double sw0;
  static double XR[4096];
  static double XI[4096];

  sw0=0;
  for (i=1;i<4096;i++) sw0+=sw[i];
  sw0/=4096.0;
  for (i=0;i<4096;i++) {
    XR[rv4096(i)]=(sw[i]-sw0)*WF2[i];
    XI[i]=0.0;
  }
  a=1;
  for (i=0;i<12;i++) {  // 4096=2^12
    for (j=0;j<4096/2/a;j++) {
      for (k=0;k<a;k++) {
        cth=cos((M_PI*k)/a);
        sth=sin((M_PI*k)/a);
        //fprintf(stderr,"%lf %lf\n",cth,sth);
        x0=XR[j*a*2+k];x1=XR[j*a*2+k+a];
        x2=XI[j*a*2+k];x3=XI[j*a*2+k+a];
        xr0=x0+x1*cth-x3*sth;
        xi0=x2+x1*sth+x3*cth;
        xr1=x0-x1*cth+x3*sth;
        xi1=x2-x1*sth-x3*cth;
        XR[j*a*2+k]=xr0;XR[j*a*2+k+a]=xr1;
        XI[j*a*2+k]=xi0;XI[j*a*2+k+a]=xi1;
      }
    }
    a*=2;
  }
  for (i=0;i<2048;i++) {
    //fprintf(stderr,"%d %lf %lf %d\n",i,XR[i],XI[i],sw[i]);
    PWR[i]=10*log10(XR[i]*XR[i]+XI[i]*XI[i]);
  }
}

int main(int argc,char *argv[]){
  int i,j,k,l,c0,c1,c2,c3,x,pw0,pw1,c,md;
  static char buf[256],path[64];
  static char s[64];
  static int data[NS];
  static double data2[NS];
  FILE *fp,*fpl;
  SpiceDouble et,spoint[3],tpoint[3],radius,lon,lat,dist,alt,dum;
  double iph,cml,io,ele;
  long int date;
  static unsigned short v[6];
  //static unsigned short v2[14];
  double f;
  static short v1[2*4096];
  static short v2[2*4096];
  static double v3[4096];
  static double v4[4096];
  static double p1[2048];
  static double p2[2048];
  double ch1;
  double ch2;
  int n1;
  int n2;
  int m1;
  int m2;
  double rr;

  sscanf(argv[1],"%ld",&date);
  sscanf(argv[3],"%d",&pw0);
  sscanf(argv[4],"%d",&pw1);
  sscanf(argv[5],"%d",&m1);
  sscanf(argv[6],"%d",&m2);
  sscanf(argv[7],"%d",&n1);
  sscanf(argv[8],"%d",&n2);
  sscanf(argv[9],"%lf",&rr);

  furnsh_c("/home/jupiter/kernel/naif0012.tls");
  furnsh_c("/home/jupiter/kernel/pck00010.tpc");
  furnsh_c("/home/jupiter/kernel/de421.bsp");
  furnsh_c("/home/jupiter/kernel/jup230l.bsp");

  fft4096_init();

  //if (w) printf("Content-type: image/png\n\n");
  KGopen(X0,Y0,"color","png",stdout);
  KGstring(16,X1,Y1,1,"Iitate Array Wideband Spectrogram",black);
  sprintf(buf,"Date:%s",argv[1]);
  KGstring(16,X2,Y1,3,buf,black);
  KGbox(X1,Y1,X2,Y2,black);
  KGbox(X1,Y2,X2,Y3,black);

  x=X1;

  for (i=0;i<1200;i++) {
    //fprintf(stderr,"%d\n",i);
    sprintf(buf,"/home/jupiter/wdp6j/1sec/%s/%s%s/%s%s%1d%02d%02d.gz",
            argv[1],argv[1],argv[2],argv[1],argv[2],i/120,(i%120)/2,50*(i%2));
    if (NULL!=(fp=fopen(buf,"r"))) {
      fclose(fp);
      sprintf(buf,"gunzip < /home/jupiter/wdp6j/1sec/%s/%s%s/%s%s%1d%02d%02d.gz",
              argv[1],argv[1],argv[2],argv[1],argv[2],i/120,(i%120)/2,50*(i%2));
      fp=popen(buf,"r");
      //fread(dat,sizeof(dat),fp);  // short(2)x900x2
      fread(v1,sizeof(v1),1,fp);
      fread(v2,sizeof(v2),1,fp);

      for (k=0;k<4096;k++) {
        if      (m1==0) v3[k]=v1[k];
        else if (m1==1) v3[k]=v2[k]*rr;
        else { 
          if (n1>=0) v3[k]=(v1[k]+v2[k+n1]*rr)/2;
          else       v3[k]=(v1[k-n1]+v2[k]*rr)/2;
        }
        if      (m2==0) v4[k]=v1[k];
        else if (m2==1) v4[k]=v2[k]*rr;
        else {
          if (n2>=0) v4[k]=(v1[k]+v2[k+n2]*rr)/2;
          else       v4[k]=(v1[k-n2]+v2[k]*rr)/2;
        }
      } 

      fft4096(v3,p1);
      fft4096(v4,p2);

      for (x=0;x<300;x++) {
        ch1=p1[452+x*3];
        ch2=p2[452+x*3];
        if (ch1<p1[452+x*3+1]) ch1=p1[452+x*3+1];
        if (ch2<p2[452+x*3+1]) ch2=p2[452+x*3+1];
        if (ch1<p1[452+x*3+2]) ch1=p1[452+x*3+2];
        if (ch2<p2[452+x*3+2]) ch2=p2[452+x*3+2];

        f=13.56+(40.56-13.56)/300.0*x;  // 452=>13.56 1352=>40.56
	f=20*log10(f);
        ch1=ch1+f-270-8;  // Ch1 -> dBW/m^2/Hz
        ch2=ch2+f-270-8;  // Ch2 -> dBW/m^2/Hz
        ch1=(ch1-pw0)/(pw1-pw0)*128;
        ch2=(ch2-pw0)/(pw1-pw0)*128;
        if (ch1<0) ch1=0;
        if (ch2<0) ch2=0;
        if (ch1>127) ch1=127;
        if (ch2>127) ch2=127;
        //fprintf(stderr,"%d %d\n",(int)ch1,(int)ch2);
	KGbox(X1+i,Y2-300-x,X1+i+1,Y2-300-x,col[(int)ch1]);
	KGbox(X1+i,Y2    -x,X1+i+1,Y2    -x,col[(int)ch2]);
      }
      pclose(fp);
    }
  }

  KGframe(X1,Y1,X2,Y2-300,black);
  KGframe(X1,Y2-300,X2,Y2,black);
  KGline(X1     ,Y3,X1     ,Y3+6,black);
  KGline(X1+ 240,Y3,X1+ 240,Y3+6,black);
  KGline(X1+ 480,Y3,X1+ 480,Y3+6,black);
  KGline(X1+ 720,Y3,X1+ 720,Y3+6,black);
  KGline(X1+ 960,Y3,X1+ 960,Y3+6,black);
  KGline(X1+1200,Y3,X1+1200,Y3+6,black);
  sprintf(s,"%c%c:%c0",argv[2][0],argv[2][1],argv[2][2]);KGstring(16,X1    ,Y3+8,8,s,black);
  sprintf(s,"%c%c:%c2",argv[2][0],argv[2][1],argv[2][2]);KGstring(16,X1+240,Y3+8,8,s,black);
  sprintf(s,"%c%c:%c4",argv[2][0],argv[2][1],argv[2][2]);KGstring(16,X1+480,Y3+8,8,s,black);
  sprintf(s,"%c%c:%c6",argv[2][0],argv[2][1],argv[2][2]);KGstring(16,X1+720,Y3+8,8,s,black);
  sprintf(s,"%c%c:%c8",argv[2][0],argv[2][1],argv[2][2]);KGstring(16,X1+960,Y3+8,8,s,black);

  KGline(X1,Y2-300-293,X1-6,Y2-300-293,black);
  KGline(X1,Y2-300-182,X1-6,Y2-300-182,black);
  KGline(X1,Y2-300- 71,X1-6,Y2-300- 71,black);
  KGline(X1,Y2-293,X1-6,Y2-293,black);
  KGline(X1,Y2-182,X1-6,Y2-182,black);
  KGline(X1,Y2- 71,X1-6,Y2- 71,black);
  KGstring(16,X1-8,Y2-300-293,6,"40",black);
  KGstring(16,X1-8,Y2-300-182,6,"30",black);
  KGstring(16,X1-8,Y2-300- 71,6,"20",black);
  KGstring(16,X1-8,Y2-293,6,"40",black);
  KGstring(16,X1-8,Y2-182,6,"30",black);
  KGstring(16,X1-8,Y2- 71,6,"20",black);

  KGstring(16,X1-36,Y1+150,-2,"RH (Array)",black);
  KGstring(16,X1-36,Y1+450,-2,"RH (Single)",black);

  for (i=0;i<1200;i++) {
    sprintf(buf,"%c%c%c%c-%c%c-%c%cT%c%c:%c%d:%02d",
	    argv[1][0],argv[1][1],argv[1][2],argv[1][3],
            argv[1][4],argv[1][5],argv[1][6],argv[1][7],
	    argv[2][0],argv[2][1],argv[2][2],i/120,(i%120)/2);
    str2et_c(buf,&et);

    subpt_c("Near point","EARTH",et,"NONE","JUPITER",spoint,&dum);
    tpoint[0]=cos(GLAT/dpr_c())*cos(GLON/dpr_c());
    tpoint[1]=cos(GLAT/dpr_c())*sin(GLON/dpr_c());
    tpoint[2]=sin(GLAT/dpr_c());
    if (spoint[0]*tpoint[0]+spoint[1]*tpoint[1]+spoint[2]*tpoint[2]>=0) 
      c=red;
    else                                                                
      c=black;
    KGline(X1+i,Y2,X1+i,Y2+10,c);

    subpt_c("Near point","JUPITER",et,"LT+S","EARTH",spoint,&dist);
    reclat_c(spoint,&radius,&lon,&lat);
    cml=-lon*dpr_c();
    while (cml<   0) cml+=360;
    while (cml>=360) cml-=360;

    subpt_c("Near point","JUPITER",et-dist/2.99792458e5,"LT+S","IO",spoint,&alt);
    reclat_c(spoint,&radius,&lon,&lat);
    iph=-lon*dpr_c();  // Io-Lon
    io=cml-iph+180;    // Io-Phase
    while (io<   0) io+=360;
    while (io>=360) io-=360;

    if      ((180<=cml)&&(cml<=300)&&(180<=io)&&(io<=260)) c=cyan;   // Io-A
    else if (( 15<=cml)&&(cml<=240)&&( 40<=io)&&(io<=110)) c=yellow; // Io-B
    else if ((200<=cml)&&(cml<=300))                       c=blue; // Non-Io-A
    else                                                   c=black;
    KGline(X1+i,Y2+10,X1+i,Y2+20,c);
  }

  KGframe(X1,Y2,X2,Y3,black);
  KGline(X1,Y2+10,X2,Y2+10,black);
  KGclose();
}