/*
 1995,1,26 own compress mothed
 2008,10 
  4096*4096*2 or *4 --> 512*512*2  --> 512*256 --> ~40000
   raw data       shrink          comp4       jcomp
                             white_black
  2880+17+~40000
*/

#include <stdio.h>
#include <math.h>
FILE *fp,*fp1;
char head[36][80];
char table[17];
char c256[131072];
char b[512*512*2];
int pa,i,j,ibit,nn,n1,n2,m1,m2,mm,k=0;
char *a,f1[80],f2[78];
unsigned char sc,spc;
char cl[12],cb[12];
float alpha,delta,sunl,sunb;

main(ac,av)
int ac; char *av[];
{
  if(ac<2){
    printf("\n\n\t *** compress fits_file ***\n");
    printf("\t\t\tde_compress command is: dh5,dh8,dh9...\n");
    printf("\n\tusage: ch554k batch_file\n\n"); exit(0);
  }
  fp1=fopen(av[1],"r");
l10:
  fgets(f2,60,fp1);
  if(feof(fp1)!=0)exit(0);
  k++; for(i=0;i<60;i++)if(f2[i]<32)f2[i]=0;	
  printf("%3d:  %s\n",k,f2);

  pa=0; sc=spc=0;
  fp=fopen(f2,"rb");
  readhead(head);
  a=(char *)malloc(nn);
  readata(a,b); 
  free(a);
  comp4(b,c256);
  jcomp(c256,b,table);   
  fclose(fp);
  j=strlen(f2);	
  for(i=j;i>0;i--)if(f2[i]=='/')break;
  strcpy(f1,"j_"); strcat(f1,&f2[i+1]);
  unlink(f1);
  fp=fopen(f1,"wb");
  strncpy(&head[1][64],"bit =  4",8);
  sprintf(&head[2][67],"%5d",pa);      head[2][72]=32;
  sprintf(&head[3][64],"m1 = %3d",m1); head[3][72]=32;
  sprintf(&head[4][64],"m2 = %3d",m2); head[4][72]=32;
  sscanf(&head[22][20],"%f",&alpha);
  sscanf(&head[23][20],"%f",&delta);
  ad_lb(alpha,delta,&sunl,&sunb);
  toms(sunl,cl,2);
  toms(sunb,cb,0);
  sprintf(&head[19][56],"l_sun = %s",cl); head[19][73]=32;
  sprintf(&head[20][56],"b_sun = %s",cb); head[20][73]=32;
  fwrite(head,1,2880,fp); fwrite(table,1,17,fp); fwrite(b,1,pa,fp);
  fclose(fp);

  goto l10;
}

readhead(head)
char head[36][80];
{
  char a[36][80];
  int i,j,k=0,ipos;
  for(i=0;i<36;i++)for(j=0;j<80;j++)head[i][j]=32;
  strncpy(&head[0][0],"SIMPLE  =",9);
  strncpy(&head[1][0],"BITPIX  =",9);
  strncpy(&head[2][0],"NAXIS   =",9);
  strncpy(&head[3][0],"NAXIS1  =",9);
  strncpy(&head[4][0],"NAXIS2  =",9);
  strncpy(&head[5][0],"CRVAL1  =",9);
  strncpy(&head[6][0],"CRVAL2  =",9);
  strncpy(&head[7][0],"DATE-OBS=",9);
  strncpy(&head[8][0],"TIME    =",9);
  strncpy(&head[9][0],"EXPOSURE=",9);
  strncpy(&head[10][0],"RA      =",9);
  strncpy(&head[11][0],"DEC     =",9);
  strncpy(&head[12][0],"HA      =",9);
  strncpy(&head[13][0],"OBJECT  =",9);
  strncpy(&head[14][0],"INSTRUME=",9);
  strncpy(&head[15][0],"OBSERVER=",9);
  strncpy(&head[16][0],"A81     =",9);
  strncpy(&head[17][0],"A82     =",9);
  strncpy(&head[18][0],"A83     =",9);
  strncpy(&head[19][0],"A84     =",9);
  strncpy(&head[20][0],"A85     =",9);
  strncpy(&head[21][0],"A86     =",9);
  strncpy(&head[22][0],"A87     =",9);
  strncpy(&head[23][0],"A88     =",9);
  strncpy(&head[24][0],"GALLONG =",9);
  strncpy(&head[25][0],"GALLATI =",9);
  strncpy(&head[26][0],"EXTINCT =",9);
  strncpy(&head[27][0],"AIRMASS =",9);
  strncpy(&head[28][0],"AZIMUTH =",9);
  strncpy(&head[29][0],"ALTITUDE=",9);
  strncpy(&head[30][0],"MPHASE  =",9);
  strncpy(&head[31][0],"MAZIMUTH=",9);
  strncpy(&head[32][0],"MALITIUD=",9);
  strncpy(&head[33][0],"MANGLE  =",9);
  strncpy(&head[34][0],"MDIRECT =",9);
  strncpy(&head[35][0],"END      ",9);
l10:
  k++;
  if(k==5){printf(" bad fits head\n"); exit(0);}
  fread(a,1,2880,fp);
  ahead(a,head,0,"SIMPLE  ");
  ahead(a,head,1,"BITPIX  ");
  ahead(a,head,2,"NAXIS   ");
  ahead(a,head,3,"NAXIS1  ");
  ahead(a,head,4,"NAXIS2  ");
  ahead(a,head,5,"CRVAL1  ");
  ahead(a,head,6,"CRVAL2  ");
  ahead(a,head,7,"DATE-OBS");
  ahead(a,head,8,"TIME    ");
  ahead(a,head,9,"EXPOSURE");
  ahead(a,head,10,"RA      ");
  ahead(a,head,11,"DEC     ");
  ahead(a,head,12,"HA      ");
  ahead(a,head,13,"OBJECT  ");
  ahead(a,head,14,"INSTRUME");
  ahead(a,head,15,"OBSERVER");
  ahead(a,head,16,"A81     ");
  ahead(a,head,17,"A82     ");
  ahead(a,head,18,"A83     ");
  ahead(a,head,19,"A84     ");
  ahead(a,head,20,"A85     ");
  ahead(a,head,21,"A86     ");
  ahead(a,head,22,"A87     ");
  ahead(a,head,23,"A88     ");
  ahead(a,head,24,"GALLONG ");
  ahead(a,head,25,"GALLATI ");
  ahead(a,head,26,"EXTINCT ");
  ahead(a,head,27,"AIRMASS ");
  ahead(a,head,28,"AZIMUTH ");
  ahead(a,head,29,"ALTITUDE");
  ahead(a,head,30,"MPHASE  ");
  ahead(a,head,31,"MAZIMUTH");
  ahead(a,head,32,"MALITIUD");
  ahead(a,head,33,"MANGLE  ");
  ahead(a,head,34,"MDIRECT ");

  ipos=indexpos(a,"NAXIS1  ",36);
  if(ipos!=36)sscanf(&a[ipos][20],"%d",&n1);

  ipos=indexpos(a,"NAXIS2  ",36);
  if(ipos!=36)sscanf(&a[ipos][20],"%d",&n2);

  ipos=indexpos(a,"BITPIX  ",36);
  if(ipos!=36){
    sscanf(&a[ipos][20],"%d",&ibit);
    nn=n1*n2*2;
    if(ibit!=16)nn=nn+nn;
  }

  ipos=indexpos(a,"IMNAME  ",36);
  if(ipos!=36){
    strncpy(&head[13][60],"IMNAME =",8); 
    strncpy(&head[13][68],&a[ipos][24],5);
  }

  ipos=indexpos(a,"VOLT1   ",36);
  if(ipos!=36){
    strncpy(&head[14][51],"TEMPERATURE =",13);
    strncpy(&head[14][64],&a[ipos][20],9);
  }

  ipos=indexpos(a,"VOLT2   ",36);
  if(ipos!=36){
    if(a[ipos][33]!='V'){
      strncpy(&head[9][60],"SKYADU =",8);
      strncpy(&head[9][68],&a[ipos][22],8);
    }
  }

  ipos=indexpos(a,"RA2000  ",36);
  if(ipos!=36){
    strncpy(&head[22][54],"RA2000  =",9);
    strncpy(&head[22][63],&a[ipos][16],13);
  }

  ipos=indexpos(a,"DEC2000 ",36);
  if(ipos!=36){
    strncpy(&head[23][54],"DEC2000 =",9);
    strncpy(&head[23][63],&a[ipos][16],13);
  }

  ipos=indexpos(a,"END     ",36);
  if(ipos==36)goto l10;

  ipos= (n1>n2) ? n1 :n2;
  mm=(ipos-1)/512+1;
  m1=n1/mm; m2=n2/mm;
}

ahead(a,head,n,c8)
char a[36][80],head[36][80],c8[8];
int n;
{
  int ipos;
  ipos=indexpos(a,c8,36);
  if(ipos!=36)strncpy(&head[n][0],&a[ipos][0],80);
}

short len[17]={3,3,3,6,6,7,7,8,8,9,9,10,10,11,11,11,11};
short dat[17]={4,5,6,0x3c,0x3d,0x7c,0x7d,0xfc,0xfd,0x1fc,0x1fd,
               0x3fc,0x3fd,0x7fc,0x7fd,0x7fe,0x7ff};
jcomp(b,a,table)
unsigned char b[],a[],table[];
{
  unsigned char chtable[256], c1[256], cc;
  int b1[256],i,j,m,n;
  short c,k;
  n=m1*m2/2;
  for(i=0;i<256;i++){c1[i]=i; b1[i]=0; chtable[i]=99; }
  for(i=0;i<n;i++)b1[b[i]]++; 
  for(i=0;i<255;i++)for(j=i+1;j<256;j++)if(b1[i]<b1[j]){
    m=b1[i]; b1[i]=b1[j]; b1[j]=m;
    cc=c1[i]; c1[i]=c1[j]; c1[j]=cc;
  }
  for(i=0;i<17;i++){ table[i]=c1[i+1]; chtable[c1[i+1]]=i; }
  for(i=0;i<n;i++){
    c=b[i];
    if(c==0)putbit(1,0,a);
    else {
      k=chtable[c];
      if(k>16){ putbit(4,14,a); putbit(8,c,a); }
      else putbit(len[k],dat[k],a); 
    }
  }
  putbit(8,0,a);       /*  adjust last byte */ 
  pa--;
}


static short d[11]={1,2,4,8,16,32,64,128,256,512,1024};
putbit(m,v,a)
char a[];
short m,v;
{
  short i;
  for(i=m-1;i>=0;i--){
    sc<<=1;
    spc++;   
    if(v&d[i])sc++;
    if(spc==8){ a[pa++]=sc; spc=sc=0; }
  }
  return;
}

readata(a,b)
char a[],b[];
{
  fread(a,1,nn,fp);
  if(ibit==16){ swap2(a,nn);   shrink2(a,b); }
  else { swap4(a,nn);  shrink4(a,b); }
}

shrink2(a,b)
short a[],b[];
{
  int i,j,k=0,i1,j1,k1;
  float xm,y;
  xm=mm*mm;
  for(i=0;i<m2*mm;i+=mm)for(j=0;j<m1*mm;j+=mm){
    y=0.;
    for(i1=0;i1<mm;i1++){
      k1=(i1+i)*n1+j;
      for(j1=0;j1<mm;j1++)y+=a[k1+j1];
    }
    b[k++]=y/xm+0.5;
  }
}

shrink4(a,b)
float a[];
short b[];
{
  int i,j,k=0,i1,j1,k1;
  float x,xm,y;
  xm=mm*mm;
  for(i=0;i<m2*mm;i+=mm)for(j=0;j<m1*mm;j+=mm){
    y=0.;
    for(i1=0;i1<mm;i1++){
      k1=(i1+i)*n1+j;
      for(j1=0;j1<mm;j1++){ x=a[k1+j1]; if(x>32767.)x=32767.; y+=x; }
    }
    b[k++]=y/xm+0.5;
  }
}


comp4(b,c)
short b[]; unsigned char c[];
{
  float white,black,d,sigma;
  int i=-1,j,j1,j2;
    
  white_(b,m1,m2,&white,&sigma);
  black=white+20.*sigma;
  d=14./(black-white);
  for(j=0;j<m1*m2/2;j++){
    i++;
    if(b[i] <= white)j1=0;
    else if(b[i] >= black)j1=14;
    else j1=(b[i]-white)*d+0.5;
    i++;
    if(b[i] <= white)j2=0;
    else if(b[i] >= black)j2=14;
    else j2=(b[i]-white)*d+0.5;
    c[j]=j1*16+j2;
  }
}


white_(map,n1,n2,peak,sigma)
short map[];
int n1,n2;
float *peak,*sigma;
{
  int *ihist;
  ihist=(int*)malloc(32767*4);
  white_1(map,n1,n2,peak,sigma,ihist);
  free(ihist);
}

white_1(map,n1,n2,xpeak,sigma,ihist)
short map[]; int n1,n2;
float *xpeak, *sigma;
long  ihist[];
{
  int i,ii;
  int imax=0;
  for(i=0;i<32767;i++)ihist[i]=0;
  for(i=0;i<n1*n2;i++){
    ii=map[i];
    ii+=101;
    if(ii<0 || ii >= 32767)continue;
    ihist[ii]++;
    if(ii>imax)imax=ii;
  }
  histt(ihist,xpeak,sigma,imax);
  *xpeak-=100.0;
}

histt(ihist,xpeak,fsigma,n)
long ihist[];
float *xpeak,*fsigma;
int   n;
{
  int ihtot=0, ihsum=0, iwid=0, maxh=0, mode=0;
  int i,ii,j,jj,k,m,nn, isigm, nfilt, noff, ixpeak, ilow, ihih,ilim,icount;
  float xpar[101],buf[201];
  float xmean=0., xmed=0, sigmsq=0., xcount=0;
  float rttpi, conv, xx, temp, xp,fs,gs,cogd,cogn,cog,xnumb,sd;

  rttpi=2.506628275;
  for(i=1;i<=n;i++){
    nn=n+1-i;
    if(ihist[nn]!=0)break;
  }
  for(i=1;i<=nn;i++)ihtot+=ihist[i];
  for(i=1;i<=nn;i++){
    ihsum+=ihist[i];
    if(ihsum < 0.1*ihtot || ihsum > 0.9*ihtot)continue;
    if(ihist[i] <= maxh)continue;
    maxh=ihist[i];
    mode=i;
  }
  for(i=1;i<=nn;i++){
    temp=ihist[i];
    xcount+=temp;
    xmean+=temp*i;
    sigmsq+=temp*i*i;
  }
  for(i=1;i<=nn;i++){
    xmed+=ihist[i];
   if(xmed > 0.5*xcount){ xmed=i; break; }
  }
  xmean/=xcount;
  gs=xcount/(maxh*rttpi);
  if(gs > 0.1*nn)gs=0.1*nn;
  isigm=gs+0.5;
  if(isigm<3)isigm=3;
  if(isigm>50)isigm=50;
  nfilt=isigm*4+1;
  noff=(nfilt+1)/2;
  conv=0.5/(float)(isigm*isigm);
  for(i=1;i<=noff;i++){
    xx=i-noff;
    buf[i]=exp(-conv*xx*xx);
    buf[nfilt+1-i]=buf[i];
  }
  xp=mode;
  for(k=1;k<=2;k++){
    ixpeak=xp+0.5;
    if(ixpeak > nn)ixpeak=nn;
    if(ixpeak <1 )ixpeak=1;
    ilow=ixpeak-isigm;
    if(ilow <1 )ilow=1;
    ihih=ixpeak+isigm;
    if(ihih >nn)ihih=nn;
    m=ihih-ilow+1;
    ii=1;  cogd=cogn=0.;
    for(i=ilow;i<=ihih;i++){
      cogd+=ihist[i];
      cogn+=ihist[i]*i;
      temp=xnumb=0.;
      for(j=1;j<=nfilt;j++){
	jj=j-noff;
	if(i+jj < 1 || i+jj > nn)continue;
	temp+=ihist[i+jj]*buf[j];
	xnumb+=buf[j];
      }
      temp = (xnumb > 0.) ? temp/=xnumb : 0. ;
      if(temp < 1.)temp=1.;
      xpar[ii++]=log((double)temp);
    }
    parbol(xpar,m,&xp,&sd);
    xx=sd*sd-isigm*isigm;
    if(xx < 0.)xx=0.;
    fs=sqrt((double)xx);
    xp+=ilow-1;
    cog=(cogd > .5) ? cogn/cogd : xp;
    xx=xp-mode; if(xx <0.)xx=-xx;
    if(xx > gs)xp=cog;
  }
  mode--;
  xp=xp-1.0;
  xmean=xmean-1.0;
  xmed=xmed-1.0;
  if(xp > nn-gs)xp=xmed;
  if(xp < gs)xp=xmed;
  ixpeak=xp+1.5;
  icount=ihist[ixpeak];
  ilim=0.68*xcount+0.5;
  for(i=1;i<=nn;i++){
    ii=ixpeak+i;
    if(ii <= nn){
      iwid++;
      icount+=ihist[ii];
      if(icount > ilim)break;
    }
    ii=ixpeak-i;
    if(ii >= 1){
      iwid++;
      icount+=ihist[ii];
      if(icount > ilim)break;
    }
  }
  gs=0.5*iwid;
  if(gs<fs)fs=gs;
  if(fs<1.)fs=1.;
  *xpeak=xp;
  *fsigma=fs;
}