#include <stdio.h>
#include <math.h>
#include <string.h>
FILE  *fp;
short bb[4064][4064];
char  head[36][80];
char  pp[80],p0[80];
float p1,p2;

short b[508*508];
char c256[508*508];
int  m1,m2;

main(ac,av)
int ac; char *av[];
{
  int ibit,n1,n2;
  if(ac<2){ printf("\n\t Usage: edo fitsname (4064*4064)\n"); exit(0); }
  fp=fopen(av[1],"rb");
  if(fp==0){printf("\n\t input file not found!\n"); exit(0); }
  fread(head,36,80,fp);
  sscanf(&head[1][22],"%d",&ibit);
  sscanf(&head[3][22],"%d",&n1);
  sscanf(&head[4][22],"%d",&n2);
  if(n1!=n2 ||n1!=4064){printf("\n\tnot 4064*4064 fits!\n"); exit(0); }
  while(index36(head,"END     ")==36)fread(head,36,80,fp);
  readata(bb,n1,ibit);
  strcpy(p0,"edo"); p1=p2=1.;  // define as r90.par
  togif();
}  

index36(head,f1)
char head[][80],f1[8];
{
  int i=-1,j;
l10:
  i++; if(i==36)return i;
  for(j=0;j<8;j++) if(head[i][j]!=f1[j])goto l10;
  return i;
}

readata(bb,n,ibit)
short bb[4064][4064];
int n,ibit;
{
  float a[4064];
  int i,j;
  if(ibit==-32){ 
    for(i=0;i<n;i++){
      fread(a,4,n,fp); swap4(a,4*n); 
      for(j=0;j<n;j++)bb[i][j]=a[j]/2;    
    }
  } else { fread(bb,2,n*n,fp); swap2(bb,2*n*n); }
  fclose(fp);
}

togif()
{
  m1=m2=508;
  shrink8(bb,b);
  comp4(b,c256);
  strcpy(pp,p0);  strcat(pp,".gif");
  printf("%s\n",pp);
  fp=fopen(pp,"wb");
  jcomp();
  fclose(fp);

  strcpy(pp,"convert "); strcat(pp,p0); strcat(pp,".gif ");
                         strcat(pp,p0); strcat(pp,".jpg");
  system(pp); printf("%s\n",pp);
  strcpy(pp,p0);  strcat(pp,".jpg");
  strcpy(pp,"rm "); strcat(pp,p0); strcat(pp,".gif");
  system(pp); printf("%s\n",pp);
// display *.jpg file
  strcpy(pp,"eog "); strcat(pp,p0); strcat(pp,".jpg &");
  system(pp); printf("%s\n",pp);
}

shrink8(a,b)
short a[4064][4064],b[508][508];
{
  int i,j,i1,j1,y;
  for(i=0;i<508;i++)for(j=0;j<508;j++){
    y=0;
    for(i1=0;i1<8;i1++)for(j1=0;j1<8;j1++)y+=a[i1+i*8][j1+j*8];
    b[i][507-j]=y/64;
  }
}

comp4(b,c)
short b[]; unsigned char c[];
{
  float white,black,d,sigma;
  short bb[508*508];
  int i,j;
  white_(b,508,508,&white,&sigma);
  black=white+3.*sigma*p2;
  white-=6.*sigma*p1;
  d=250./(black-white);
  for(j=0;j<508*508;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;
int  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)
int 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 int TempCode,jj;
jcomp()
{	
  int i;
  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=(m1-Rows)*m1;
    for(i=0;i<m1;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()
{
  int i;
  Encode = 258;  RunBits = 9;  MaxCodeSize =512;
  for(i = 0; i < TABLESIZES; i++) EncodeTable[i] = 0;
}

FillBlockBuf()
{
  TempCode |= (unsigned int)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;
}

swap2(char *a, int n)
{
  register char *p;
  register int i;
  register char c;
  p = a;
  for(i=0;i<n;i+=2){
    c = *p;
    *p++ = *(p+1);
    *p++ = c;
  }
}

swap4(char *a, int n)       /*     int a[n] or float a[n]    */
{
  register char *p;
  register int i;
  register char c,d;
  p = a;
  for(i=0;i<n;i+=4){
    c = *p;
    *p++ = *(p+3);
    d = *p;
    *p++ = *(p+1);
    *p++ = d;
    *p++ = c;
  }
}

parbol(a,n,xpeak,sd)
float a[],*xpeak,*sd; int n;
{
  float aa,fi,fi2;
  float cmatx[25][25],bvect[25];   /* same as fortran , 1--24 */
  float suma=0.,sumay=0.,sumayy=0.,sumy1=0.,sumy2=0.,sumy3=0.,sumy4=0.;
  int i,j;

  for(i=1;i<25;i++){
    bvect[i]=0.;
    for(j=1;j<25;j++)cmatx[i][j]=0.;
  }
  for(i=1;i<=n;i++){
    aa=a[i]; fi=i; fi2=i*i;
    suma+=aa;  sumay+=aa*fi; sumayy+=aa*fi2;
    sumy1+=fi; sumy2+=fi2;   sumy3+=fi*fi2;   sumy4+=fi2*fi2;
  }
  cmatx[1][1]=n;     cmatx[1][2]=sumy1; cmatx[1][3]=sumy2;
  cmatx[2][1]=sumy1; cmatx[2][2]=sumy2; cmatx[2][3]=sumy3;
  cmatx[3][1]=sumy2; cmatx[3][2]=sumy3; cmatx[3][3]=sumy4;
  bvect[1]=suma;     bvect[2]=sumay;    bvect[3]=sumayy;
  solve(cmatx,bvect,3);
  if(bvect[3] == 0) *xpeak= *sd=0.;
  else {
    aa=bvect[3]*2.;
    *xpeak=-bvect[2]/aa;
    if(aa < 0.)aa=-aa;
    *sd=sqrt((double)1.0/aa);
  }
}

solve(a,b,m)
float a[25][25],b[25]; int m;
{
  int i,j,k,l;
  float big,temp;

  for(i=1;i<m;i++){
    big=0.;
    for(k=i;k<=m;k++){
      temp=a[k][i];
      if(temp < 0.)temp=-temp;
      if(temp > big){ big=temp; l=k; }
    }
    if(big == 0.){
      for(k=1;i<=m;k++)b[k]=0;
      printf(" Zero determinant from SOLVE\n");
      return;
    }
    if(i != l){
      for(j=1;j<=m;j++){ temp=a[i][j]; a[i][j]=a[l][j]; a[l][j]=temp;}
      temp=b[i]; b[i]=b[l]; b[l]=temp;
    }
    big=a[i][i];
    for(j=i+1;j<=m;j++){
      temp=a[j][i]/big;
      b[j]-=temp*b[i];
      for(k=i;k<=m;k++)a[j][k]-=temp*a[i][k];
    }
  }
  for(i=1;i<=m;i++){
    l=m+1-i;
    if(a[l][l]==0.){ b[l]=0.; continue; }
    temp=b[l];
    if(l != m)for(j=2;j<=i;j++){ k=m+2-j; temp-=a[l][k]*b[k]; }
    b[l]=temp/a[l][l];
  }
}