/*
 1995,1,26 own compress mothed
  2048*2048*2 or *4 --> 512*512*2  --> 512*512 --> 
   raw data       shrink          comp4       gif
                             white_black
*/
#include <stdio.h>
#include <math.h>
FILE *fp;
char head[36][80];
char c256[512*512];
char b[512*512*2];
int  i,ibit,nn,n1,n2,m1,m2,mm,k=0;
//unsigned char sc,spc;
char cl[12],cb[12];
float alpha,delta,sunl,sunb;

main(ac,av)
int ac; char *av[];
{
  char *a;
  char f1[60];
  if(ac<3){
    printf("\n\n\t *** change fits_file to GIF (512*512*256) ***\n");
    printf("\n\t Ex:   gh6 in(p*.fit) out(*.gif)  2004,8\n\n");
    exit(0);
  }
  fp=fopen(av[1],"rb");
  readhead(head);
  a=(char *)malloc(nn);
  readata(a,b); 
  free(a);
  comp4(b,c256);
  fclose(fp);
  fp=fopen(av[2],"wb");
  strncpy(&head[1][64],"bit =  8",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;
  jcomp();
  fclose(fp);
  strcpy(f1,av[2]);
  for(i=0;i<strlen(f1);i++)
	  if(f1[i]=='.')break;
  f1[i]=0;
  strcat(f1,".txt");
  fopen(f1,"w");
  for(i=0;i<36;i++){
    head[i][79]=0;
    fprintf(fp,"%s\n",&head[i][0]);
  }
  fclose(fp);
}

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);
}

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,j;
  white_(b,m1,m2,&white,&sigma);
  black=white+20.*sigma;
  d=250./(black-white);
  for(j=0;j<m1*m2;j++){
    i=(b[j]-white)*d+0.5;
    if(i<0)i=0;
    if(i>255)i=255;
    c[j]=i;
  }
}

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;
}

#define MAXCODES 4096
#define TABLESIZES 4999
unsigned char SuffixTable[4999], ByteBuf[2592], BlockBuf[256];
int Encode, RunBits, MaxCodeSize, ByteCount, ShiftBits;
int Code, Width, Height, Dots, Rows;
int PrefixCode, SuffixCode, RandomIndex, Val;
int EncodeTable[4999], PrefixTable[4999];
short head1[7]={0x4947,0x3846,0x6137,0,0,0xf7,0};
short head2[5]={0,0,0,0,0x800};
unsigned long TempCode,jj;
jcomp()
{	
  head1[3]=head2[2]=Width=m1;
  head1[4]=head2[3]=Height=m2;
  fwrite(head1,1,13,fp);
  for(i=0;i<256;i++){
    ByteBuf[0]=ByteBuf[1]=ByteBuf[2]=i;
    fwrite(ByteBuf,1,3,fp);
  }
  fputc(0x2c,fp);
  fwrite(head2,2,5,fp);
  ByteCount = ShiftBits = TempCode = 0;
  ClearTable();
  Code = 256;
  FillBlockBuf();
  for(Rows = 0; Rows < Height; Rows++)  {
	  jj=(512-Rows)*512;
    for(i=0;i<512;i++)ByteBuf[i]=c256[--jj];	  
    if(Rows == 0) { PrefixCode = ByteBuf[0];  Dots = 1; }
    else Dots = 0;
    while( Dots < Width)  {
      SuffixCode = ByteBuf[Dots++];
      RandomIndex = PrefixCode ^ (SuffixCode << 4);
      if(RandomIndex == 0)  Val = 1;
      else Val = TABLESIZES - RandomIndex;
      while(1)  {
        if(EncodeTable[RandomIndex] == 0)  {
          Code = PrefixCode;
          FillBlockBuf();
          if(Encode == MAXCODES)  {
            Code = 256;
            FillBlockBuf();
            ClearTable();
          }
          else  {
            if(Encode == MaxCodeSize)  {
              MaxCodeSize <<= 1;
              RunBits++;
            }
            PrefixTable[RandomIndex] = PrefixCode;
            SuffixTable[RandomIndex] = SuffixCode;
            EncodeTable[RandomIndex] = Encode++;
          }
          PrefixCode = SuffixCode;
          break;
        }
        if(PrefixTable[RandomIndex] == PrefixCode &&
           SuffixTable[RandomIndex] == SuffixCode)  {
          PrefixCode = EncodeTable[RandomIndex];
          break;
        }
        else  {
          RandomIndex -= Val;
          if(RandomIndex < 0)  RandomIndex += TABLESIZES;
        }
      }
    }
  }
  Code = PrefixCode;
  FillBlockBuf();
  Code = 257;
  FillBlockBuf();
  if(ShiftBits > 0 || ByteCount > 0)   {
    BlockBuf[++ByteCount] = TempCode & 0x00FF;
    SaveCode();
  }
  fputc(0, fp);  fputc(';', fp);
}

ClearTable()
{
  Encode = 258;  RunBits = 9;  MaxCodeSize =512;
  for(i = 0; i < TABLESIZES; i++) EncodeTable[i] = 0;
}

FillBlockBuf()
{
  TempCode |= (unsigned long)Code << ShiftBits;
  ShiftBits += RunBits;
  while(ShiftBits >= 8)  {
    BlockBuf[++ByteCount] = TempCode & 0x00FF;
    if(ByteCount == 255)  SaveCode();
    TempCode >>= 8;
    ShiftBits -= 8;
  }
}

SaveCode()
{
  BlockBuf[0] = ByteCount;
  fwrite(BlockBuf, 1, ByteCount+1, fp);
  ByteCount = 0;
}