#include <stdio.h>
#include <math.h>
// 2003,4,15
// 2003,12,23
// 2005,11,22
// 2006,2,23
// 2006,3
// 2006,12
main()
{
  float ra,de;
  char rc[12],dc[12],w1[30],a[99],head[72][80],c;
  char aa[120];
  char head2[72][80];
  char color,observer[30],ccd_num[12];
  FILE *fp,*fp1,*fp0;
  int run_number,i,j,k,w3,w2[30],im2,im3;
  int iy,im,id,th,tm,jd;
  float ts,epoch,epochold;
  int mon[12]={31,28,31,30,31,30,31,31,30,31,30,31};
  int i1,i2,j1,j2,j3;
  float x,y;
  int n1,n2,c1;
  char mm[80],m4[9];
  char f1[20],obsfile[24],ff[20];
  char f2[20],ach[50];
  short oneline[4160];

  puthead2(head2);                        /* init our fits_head */


  fp=fopen("/home/ccdev/RUN_NUMBER","r");
  if(fp==0){ printf("\n\tERROR: cannot found file: RUN_NUMBER\n"); exit(0); }
  fscanf(fp,"%d",&run_number); while(fgetc(fp)!=10);
  fscanf(fp,"%s",&ccd_num);
  fclose(fp);

  fp=fopen("/home/ccdev/FILTER","r");
  if(fp==0){ printf("\n\tERROR: cannot found file: FILTER\n"); exit(0); }
  i=0;
l10:
  fscanf(fp,"%c %d\n",&w1[i],&w2[i]); if(feof(fp)!=0)goto l11;
  i++; if(i<30)goto l10;
l11:
  w3=i+1;
  fclose(fp);

  fp0=fopen("/home/ccdev/jiang/exp.dat","r");
/* ************** read origion fits head */
  fgets(ff,20,fp0);
  ff[strlen(ff)-1]=0;
  strcpy(f1,"/data/d");
  strcat(f1,ff);
  strcat(f1,".ccd"); 
  printf("%s:  ",f1);
  fp=fopen(f1,"rb");
  if(fp==0){ printf("********ERROR: cannot found file !\n",f1); exit(0); }
  fread(head,72,80,fp);
//  i=indexpos(head,"END     ",36);
//  if(i!=22){ printf("********ERROR: not a raw ccd image file !\n"); exit(0);}
    strcpy(f2,"/data/D");
    strcat(f2,ff);
    strcat(f2,".ccd"); 
    strcpy(ach,"mv "); strcat(ach,f1); strcat(ach," "); strcat(ach,f2);
    system(ach);

  strncpy(&head2[26][11],ccd_num,strlen(ccd_num));           // put ccd_number
  strncpy(&head2[30][11],f1,strlen(f1));               // put *.ccd name 

/* ********* read color */
  fgets(ff,20,fp0);
  sscanf(&ff[6],"%d",&j);
  color='x';
  for(i=0;i<w3;i++)if(j==w2[i])color=w1[i];
  if(color!='x'){
    for(i=0;i<w3;i++)if(color==w1[i])goto l15; 
    printf("********ERROR: cannot find color in FILTER file !\n"); exit(0);
l15:
    sprintf(&head2[30][24],"%5d%c",w2[i],color); head2[30][30]=32;  // put wave
  }

/* ********* read name */
  fgets(ff,20,fp0);
  im3=0;
  if(strcmp(ff,"FLAT")==0)im3=-3;
  if(strcmp(ff,"DARK")==0){ im3=-2; color='x'; }
  if(strcmp(ff,"BIAS")==0){ im3=-1; color='x'; }
  for(i=0;i<80;i++)mm[i]=32;
  strcpy(m4,ff); m4[4]=0;
  strcpy(mm,m4);
  if(im3==-3)strncpy(&mm[32],"FLAT",4);
  if(im3==-2)strncpy(&mm[32],"DARK",4);
  if(im3==-1)strncpy(&mm[32],"BIAS",4);

/* ************* read ra dec object */
  
  fgets(ff,20,fp0);
    strcpy(rc,ff); rerang(rc,1); strcat(mm," "); strcat(mm,rc);
  fgets(ff,20,fp0);
    strcpy(dc,ff); rerang(dc,0); strcat(mm," "); strcat(mm,dc);
  fgets(ff,20,fp0);
    for(i=36;i<80;i++)mm[i]=32; strcpy(&mm[36],ff); 
  fclose(fp0);

/* *************read av[] end */

  i=indexpos(head,"DATE-END",36);           
  if(i==36){ printf(" keywords DATE-END not found!\n"); exit(0); }
  strncpy(a,&head[i][10],23);
  for(i=0;i<=23;i++)if(a[i]<'.' || a[i]==':' || a[i]>'@')a[i]=32;
  sscanf(a,"%d %d %d %d %d %f",&iy,&im,&id,&th,&tm,&ts);
// if th<8 th-=8; id-=1, redo_date, do jd
  i=indexpos(head,"EXPTIME ",36);           
  sscanf(&head[i][16],"%f",&x);                 // get exptime
  y=(th-8)*3600.+tm*60.+ts-x;
  k=0; if(y<0.){ y+=(24*3600.); k=1; } 
  i=y+0.05;                              // for exchange error
  ts=y-i;                                // for .second     
  th=i/3600; i-=(th*3600);
  tm=i/60;   i-=(tm*60);
  ts+=i; 
  if(k==1){
    id-=1;
    if(id==0){ 
      im-=1;
      if(im==0){
        iy-=1;
        im=12;
      }
      id=mon[im-1];
      if(iy%4==0 && im==2)id++;
    }
  }    
  jd=dojd(iy,im,id,th);
//  printf("%d %d %d %d %d\n",iy,im,id,th,jd);

// put obs???
  sprintf(obsfile,"/home/ccdev/obs%04d",jd); strcat(obsfile,".txt");
//  printf("%s\n",obsfile);
  
  fp1=fopen(obsfile,"r");
  if(fp1==0){ printf("\n\tERROR: cannot found: %s\n",obsfile); exit(0); }
  fgets(observer,30,fp1); fclose(fp1);
  strncpy(&head2[31][11],observer,strlen(observer)-1);       // put observer


  c=0x2f;    
  sprintf(&head2[11][21],"%02d%c%02d%c%02d",id,c,im,c,iy-2000); 
           head2[11][10]=0x27;  head2[11][29]=0x27;               // put date
  c=':';            
  sprintf(&head2[12][19],"%02d%c%02d%c%04.1f",th,c,tm,c,ts);
           head2[12][10]=0x27;  head2[12][29]=0x27;               // put UT
  j=(iy+((im-1)*30.4+id)/365.)*10.+0.5;  epoch=j*0.1;                   
  sprintf(&head2[17][23],"%6.1f",epoch);
           head2[17][10]=0x27;  head2[17][29]=0x27;               // put epoch
  if(im3>-4){                                        // from NAME
    c=mm[17];
    mm[7]=mm[10]=mm[17]=mm[20]=mm[23]=32;
    epochold=epoch;
    sscanf(&mm[5],"%d%d%f%d%d%f%f",&i1,&i2,&x,&j1,&j2,&y,&epochold);
    x=i1+i2/60.+x/3600.;
    y=j1+j2/60.+y/3600.;
    if(c=='-')y=-y;
    astprs(x,y,epochold,&ra,&de,epoch);
    if(im3<0)ra=de=0.;
    toms1(ra,rc,1); toms1(de,dc,0);
  }

  strncpy(&head2[14][18],rc,11);              // put RA
  strncpy(&head2[15][18],dc,11);              // put DEC
  for(i=32;i<80;i++)if(mm[i]<32)mm[i]=32;
  strncpy(&head2[25][11],&mm[36],42);         // put OBJ

// put head-->head2
  sscanf(&head[3][24],"%d",&n1); 
  sscanf(&head[4][24],"%d",&n2); 
  sscanf(&head[5][24],"%d",&c1); 
  strncpy(&head2[4][0],&head[4][0],80);       // put NAXIS2
  sprintf(&head2[3][24],"%6d",n2);
  head2[3][30]=32;                            // put NAXIS1 as NAXIS2
  sprintf(&head2[30][33],"32*2 column overscan");
  head2[30][53]=32;   head2[30][74]=0x27;     // put overscan
//    if(n2==512){
//      head2[8][29]='4';
//      head2[9][29]='4';
//    }
  strncpy(&head2[10][25],m4,4);            // put IMNAME
  sscanf(&head[10][16],"%f",&ts);
  k=ts+0.5;
  sprintf(&head2[13][24],"%6d",k);
  head2[13][30]=32;                           // put exposure time
  sprintf(&head2[43][22],"%8.2f",ts);
  head2[43][30]=32;                           // put exposure time

// 20031223
  strncpy(&head2[ 6][0],&head[ 5][0],80);     // put sc
  strncpy(&head2[ 7][0],&head[ 6][0],80);     // put sr
  strncpy(&head2[18][15],&head[42][15],65);   // put temperature
  strncpy(&head2[21][0],&head[52][0],80);     // put MPP
//  strncpy(&head2[22][0],&head[14][0],80);     // put SHUTTER
  strncpy(&head2[27][0],&head[55][0],80);     // put READ-SPD
  strncpy(&head2[29][0],&head[58][0],80);     // put GAIN

  sscanf(&f1[7],"%d",&i);
  sprintf(a,"/data/d%03d%04d%s%c%03d.fit",run_number,jd,m4,color,i);
  printf("   --> %s\n",a);

  system("date >> /u/ccdev/smtlib/obs.log");
  fp1=fopen("n0.tmp","w");
  fprintf(fp1,"%s   ---> %s\n",f1,a);
  fclose(fp1);
  system("cat n0.tmp >> /u/ccdev/smtlib/obs.log");

  fp1=fopen(a,"wb");
  fwrite(head2,72,80,fp1);
  for(i=0;i<n2;i++){
    fread(oneline,2,n1,fp);
    if((n1-n2)==32)baseline1(oneline,n1,n2);
    else           baseline(oneline,n1,n2,c1);
    fwrite(oneline,2,n2,fp1);
  }
  fclose(fp);
  fclose(fp1);
  strcpy(aa,"/u/ccdev/jiang/gh7 "); strcat(aa,a);
  printf("%s\n",aa);
  system(aa);
}

baseline(a,n1,n2,c1)
short a[];
int n1,n2,c1;
{
  int i,base1,base2;
  swap2(a,n1*2);
  base1=median(&a[n1-64],32)+32000;
  base2=median(&a[n1-32],32)+32000;
  for(i=0;i<n2;i++)if((i+c1)<2048) a[i]-=base1; else a[i]-=base2;
  i=2*32000-32768;
  if(c1<2048)a[0]=base1-i;                  // left overscan writedown
  if((c1+n2)>2048)a[n2-1]=base2-i;          // right overscan writedown
//  printf("%d %d\n",base1,base2);
  swap2(a,n2*2);
} 

baseline1(a,n1,n2)
short a[];
int n1,n2;
{
  int i,base;
  swap2(a,n1*2);
  base=median(&a[n1-32],32)+32000;
  for(i=0;i<n2;i++) a[i]-=base;
  i=2*32000-32768;
  a[n2-1]=base-i;          // right overscan writedown 
  swap2(a,n2*2);
}

median(a,n)
short a[];
int   n;
{
  int i,j;
  short k;
  for(i=0;i<n-1;i++)for(j=i+1;j<n;j++)if(a[i]>a[j]){
    k=a[j]; a[j]=a[i]; a[i]=k;
  }
  i=a[n/2];
  return(i);
}


puthead2(a)
char a[72][80];
{
  int i,j;
  for(i=0;i<72;i++)for(j=0;j<80;j++)a[i][j]=32;

  strncpy(&a[ 0][0],"SIMPLE  =                    T / NORMAL FITS IMAGE            ",62);
  strncpy(&a[ 1][0],"BITPIX  =                   16 / DATA PRECISION               ",62);
  strncpy(&a[ 2][0],"NAXIS   =                    2 / NUMBER OF IMAGE DIMENSIONS   ",62);
  strncpy(&a[ 3][0],"NAXIS1  =                      / NUMBER OF COLUMNS            ",62);
  strncpy(&a[ 4][0],"NAXIS2  =                      / NUMBER OF ROWS               ",62);
  strncpy(&a[ 5][0],"BLOCKED =                    T / CHECK FOR POSSIBLE BLOCKING  ",62);
  strncpy(&a[ 6][0],"CRVAL1  =                    0 / COLUMN ORIGIN                ",62);
  strncpy(&a[ 7][0],"CRVAL2  =                    0 / ROW ORIGIN                   ",62);
  strncpy(&a[ 8][0],"CDELT1  =                    1 / COLUMN BINNING               ",62);
  strncpy(&a[ 9][0],"CDELT2  =                    1 / ROW BINNING                  ",62);
  strncpy(&a[10][0],"IMNAME  = '                  ' / OBSERVATION SEQUENCE NAME    ",62);
  strncpy(&a[11][0],"DATE-OBS= '                  ' / DATE OF START OF OBSERVATION ",62);
  strncpy(&a[12][0],"TIME    = '                  ' / UT TIME                      ",62);
  strncpy(&a[13][0],"EXPOSURE=                      / EXPOSURE TIME (SEC)          ",62);
  strncpy(&a[14][0],"RA      = '                  ' / RIGHT ASCENSION              ",62);
  strncpy(&a[15][0],"DEC     = '                  ' / DECLINATION                  ",62);
  strncpy(&a[16][0],"HA      = '                  ' / HOUR ANGLE                   ",62);
  strncpy(&a[17][0],"EPOCH   = '                  ' / EPOCH OF RA AND DEC          ",62);
  strncpy(&a[18][0],"VOLT1   =                      / TEMPERATURE                  ",62);
  strncpy(&a[19][0],"VOLT2   =                      /                              ",62);
  strncpy(&a[20][0],"VOLT3   =                      /                              ",62);
  strncpy(&a[21][0],"VOLT4   =                      /                              ",62);
  strncpy(&a[22][0],"VOLT5   =                      /                              ",62);
//21  ...read_spo
//22  ...tempcon
  strncpy(&a[23][0],"VOLT6   =                      /                              ",62);
  strncpy(&a[24][0],"VOLT7   =                      /                              ",62);
  strncpy(&a[25][0],"OBJECT  = '",11);                                    a[25][74]=0x27;
  strncpy(&a[26][0],"INSTRUME= '",11);                                    a[26][74]=0x27;
//27  ...Mpp 
  strncpy(&a[28][0],"STATUS  = 'RAW               ' / IMAGE PROCESSING STATE       ",62);
  strncpy(&a[29][0],"SHUTTER = 'OPEN              ' / SHUTTER STATE                ",62);
  strncpy(&a[30][0],"RAWMESS = '",11);                                    a[25][74]=0x27;
  strncpy(&a[31][0],"OBSERVER= '",11);                                    a[26][74]=0x27;
  strncpy(&a[32][0],"BZERO   =             32000.00 / ZERO POINT BIAS              ",62);
  strncpy(&a[33][0],"BSCALE  =             1.000000 / IMAGE SCALE FACTOR           ",62);
  strncpy(&a[34][0],"LICK    = 'FITS2             ' / SPECIAL LICK FLAG            ",62);
  strncpy(&a[35][0],"CRPIX1  =                    1 / PIXEL CORRESPONDING TO CRVAL1",62); 
  strncpy(&a[36][0],"CRPIX2  =                    1 / PIXEL CORRESPONDING TO CRVAL2",62);
  strncpy(&a[37][0],"COMMENT = '* BAO Schmidt CCD Data-Taking system",47);a[37][74]=0x27;
  strncpy(&a[38][0],"COMMENT = '",11);                                    a[38][74]=0x27;
  strncpy(&a[39][0],"COMMENT = '",11);                                    a[39][74]=0x27;
  strncpy(&a[40][0],"COMMENT = '",11);                                    a[40][74]=0x27;
  strncpy(&a[41][0],"COMMENT = '",11);                                    a[41][74]=0x27;
  strncpy(&a[42][0],"TELESCOP= 'UNKNOWN           ' / TELESCOPE USED               ",62);
  strncpy(&a[43][0],"TTIME   =                      / REAL EXPOSURE TIME IN FLOAT  ",62);
  strncpy(&a[71][0],"END                            / END OF HEADER                ",62);
}

rerang(a,k)
char a[]; int k;
{
  char c=0;
  int i,i1,i2;
  float x;
  
  for(i=0;i<12;i++){
    if(a[i]==0)for(i1=i+1;i1<12;i1++)a[i1]=32;
    if(a[i]=='-')c=1;
    if(a[i]<'.')a[i]=32;
    if(a[i]==':')a[i]=32;
  }
  i1=i2=0; x=0.;
  sscanf(a,"%d%d%f",&i1,&i2,&x);
  if(k==1)sprintf(a,"%02d:%02d:%05.2f\0",i1,i2,x);
  if(k==0)sprintf(a," %02d:%02d:%04.1f\0",i1,i2,x);
  if(c==1)a[0]='-';
}

dojd(iy,im,id,th)
int iy,im,id,th;
{
  int i,j;
  i=(im+9)*0.09;
  i=(i+iy-1900)*1.75+0.01;
  j=30.56*im;
  j=(iy-1950)*367+id+j-i+3339;
// define bj_time 18:00 new day begin
  if((th+8)<18)j--;
  j%=10000;
  return(j);
}

astprs (ra1, dec1, ep1, ra2, dec2, ep2)
/* alpha in hour, decl.in degree*/
float ra1,dec1,ep1,*ra2,*dec2,ep2;
{
  double epoch1,epoch2,r0[3],r1[3],p[3][3],arc;
  float r2,d2;
  epoch1=ep1; epoch2=ep2;
  arc=45./atan(1.);
  *ra2=ra1; *dec2=dec1;
  if(epoch1 == epoch2)return;
  r2 =ra1*15./arc; d2 =dec1/arc;
  r0[0]=cos(r2)*cos(d2); r0[1]=sin(r2)*cos(d2); r0[2]=sin(d2);
  if(epoch1 != 2000){
   astrox(epoch1, p);
   r1[0] = p[0][0] * r0[0] + p[0][1] * r0[1] + p[0][2] * r0[2];
   r1[1] = p[1][0] * r0[0] + p[1][1] * r0[1] + p[1][2] * r0[2];
   r1[2] = p[2][0] * r0[0] + p[2][1] * r0[1] + p[2][2] * r0[2];
   r0[0] = r1[0]; r0[1] = r1[1]; r0[2] = r1[2];
  }
  if(epoch2 != 2000){
    astrox(epoch2, p);
    r1[0] = p[0][0] * r0[0] + p[1][0] * r0[1] + p[2][0] * r0[2];
    r1[1] = p[0][1] * r0[0] + p[1][1] * r0[1] + p[2][1] * r0[2];
    r1[2] = p[0][2] * r0[0] + p[1][2] * r0[1] + p[2][2] * r0[2];
    r0[0] = r1[0];    r0[1] = r1[1];    r0[2] = r1[2];
  }
  *ra2  = atan2(r0[1], r0[0])/15.*arc;
  *dec2 = asin(r0[2])*arc;
  if(*ra2<0)*ra2+=24;
}

astrox (epoch, p)
double epoch,p[3][3];
{
  double t,a,b,c,ca,cb,cc,sa,sb,sc,arc;
  arc=45./atan(1.);
  astjuy(epoch,&t);
  t = (t - 2451545.0) / 36525.;
  a = t * (0.6406161 + t * (0.0000839 + t * 0.0000050));
  b = t * (0.6406161 + t * (0.0003041 + t * 0.0000051));
  c = t * (0.5567530 - t * (0.0001185 + t * 0.0000116));
  ca = cos (a/arc);
  sa = sin (a/arc);
  cb = cos (b/arc);
  sb = sin (b/arc);
  cc = cos (c/arc);
  sc = sin (c/arc);
  p[0][0] = ca * cb * cc - sa * sb;
  p[1][0] = -sa * cb * cc - ca * sb;
  p[2][0] = -cb * sc;
  p[0][1] = ca * sb * cc + sa * cb;
  p[1][1] = -sa * sb * cc + ca * cb;
  p[2][1] = -sb * sc;
  p[0][2] = ca * sc;
  p[1][2] = -sa * sc;
  p[2][2] = cc;
}

astjuy (epoch,t)
double epoch,*t;
{
  double jd;
  int year,centuy;
  year = epoch - 1;
  centuy = year / 100;
  jd =1721425.5+365.*year-centuy+ year/4 + centuy/4;
  year=epoch;
  *t= jd + (epoch - year) * 365.25;
}

toms1(aa,cc,k)
/*
c hour or degree to char_line
c if k=1 (hour case) in **:**:**.*
c    k=0 (degree       -**:**:**
c    k=2 (degree       ***:**:**
*/
char cc[12];
float aa;
int k;
{
  int i,j,m,n,l;
  float a,x;
  a=aa; if(a < 0.)a=-a;
  i=a; x=(a-i)*60.;
  j=x; x=(x-j)*60.;
  m=x;
  if(k!=0){
  n=(x-m)*100.+0.05;
  if(n >= 100) { m++;  n=0; }
  } else {
  n=(x-m)*10.+0.5;
  if(n >= 10) { m++;  n=0; }
  }
  if(m == 60) { j++;  m=0; }
  if(j == 60) { i++;  j=0; }
  if(k ==2)sprintf(cc,"%3d:%02d:%02d.%02d\0",i,j,m,n);
  if(k ==1)sprintf(cc,"%02d:%02d:%02d.%02d\0",i,j,m,n);
  if(k ==0)sprintf(cc," %02d:%02d:%02d.%d\0",i,j,m,n);
  if(k ==0)cc[0]=(aa<0.)?'-':' ';
}