#include <stdio.h>
#include <math.h>
#define scale  0.5           // extract from *.fit in step of 10 aecsec/pixel
#define maxpix 500000000
int n1,n2;
char   head[72][80],b[144],f1[30];
float  a[maxpix];      
float  c[4032][4096];
float  xx[5],yy[5];
char   c1[14],c2[14];
double a8[8],b8[8],s8[8],t8[8],aa[4],dd[4];
double rc,dc,hs2, pi,cx,cy;
int    ix[4],iy[4];
int  ccdno;  float white;
FILE   *fp,*fp1;
int    q1d[4];                      // ={ 1,12,7,10 };
int    q2u[4];                      // ={ 10,38,1,0 };
main(int ac, char **av)
{
  int    i,j,i1,i2,j1,j2;
  double x,y,alpha,delta,x1,y1;
  int    ix0,ix1,iy0,iy1;
  float  z0,z1,z2,z3;
  float  zx,zy,wx,wy;
  float  q1,q2;
  char   c9=34;

  if(ac<5){
    printf("\n\t ***** REarrange UBAND mosaic fits (add only) ********");
    printf("\n\t                                  jiang 2011,2,9\n");
    printf("\n\tUsage:  z_6 outfile RA DEC X Y");
    printf("\n\t   or:  z_6 outfile RA DEC R\n");
    printf("\n\tSample: z_6 aa.fit 23:23:23.4 +01:12:00 3200 2400");
    printf("\n\t   or:  z_6 ab.fit 02:22 2:00 120\n");
    printf("\n\t    RA DEC: (2000.0) center of pic of outfile");
    printf("\n\t    X,Y in pixel;  R (radius) in arcmin\n");
    printf("\n\tnotice: Scale fixed to 0.%c5 /pixel.",c9);
    printf("\n\t        can erase bottom & top of each ccd's bad lines");
    printf("\n\t        copy /vega/rhbin/ubandlines.dat to default'disk\n\n");
    exit(0);
  }

  fp=fopen("ubandlines.dat","r");
  if(fp){
    for(i=0;i<4;i++){ fgets(b,80,fp); sscanf(b,"%d %d",&q1d[i],&q2u[i]);}
    printf("ubandlines.dat: ");
    for(i=0;i<4;i++)printf("%d %d, ",q1d[i],q2u[i]);  printf("\n");
    fclose(fp); }
  else printf("\nNot found 'ubandlines.dat',copy it from /vega2/rhbin\n");
 
  strcpy(f1,av[1]);                     // output file name
  strcpy(c1,av[2]);  chs(c1,&rc);       // center of RA
  strcpy(c2,av[3]);  chs(c2,&dc);       //        of Dec
  sscanf(av[4],"%d",&i);  j=0;
  if(ac>5)sscanf(av[5],"%d",&j);
  if(j){ n1=i; n2=j; }
  else   n1=n2=i/scale*60.+0.5;
  if(n1*n2>maxpix){ printf("\n error: array is too large ! %d*%d\n",n1,n2);
                    exit(0); }
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  for(j=0;j<n2*n1;j++)a[j]=0.;
  toms2(rc,c1,1); toms2(dc,c2,0);
  printf("center:  %s %s       %d*%d pixels\n",c1,c2,n1,n2);
  rc*=cx; dc*=cy; a8[6]=rc; a8[7]=dc;  
  hs2=(3600.*180./pi)/scale;                 
/*
  a8[0]=-1/hs2; a8[1]=a8[2]=0; a8[3]=-a8[0]; a8[4]=n2/2./hs2; a8[5]=-n2/2./hs2;
  for(i=0;i<8;i++)printf("%14.8le\n",a8[i]); printf("\n");
*/
  b8[0]=-hs2; b8[1]=b8[2]=0; b8[3]=-b8[0]; b8[4]=n1/2; b8[5]=n2/2;
  xytoad(b8,a8);
//  for(i=0;i<8;i++)printf("%14.8le\n",a8[i]); printf("\n");
  x=n1/2.; y=n2/2.;  xy_rad(rc,dc,x,y,&alpha,&delta,1,0);
  toms2(alpha,c1,1); toms2(delta,c2,0);
  rad_xy(rc,dc,alpha,delta,&x1,&y1,0,0);
  printf("%8.1lf %8.1lf-> %s %s-> %8.1lf %8.1lf\n",x,y,c1,c2,x1,y1);
  printf("4 corner's coodination:\n");
  x=n1;    y=n2;     xy_rad(rc,dc,x,y,&alpha,&delta,1,0);
  toms2(alpha,c1,1); toms2(delta,c2,0); 
  rad_xy(rc,dc,alpha,delta,&x1,&y1,1,0);
  printf("%8.1lf %8.1lf-> %s %s-> %8.1lf %8.1lf\n",x,y,c1,c2,x1,y1);

  x=0;     y=n2;     xy_rad(rc,dc,x,y,&alpha,&delta,1,0);
  toms2(alpha,c1,1); toms2(delta,c2,0);
  rad_xy(rc,dc,alpha,delta,&x1,&y1,1,0);
  printf("%8.1lf %8.1lf-> %s %s-> %8.1lf %8.1lf\n",x,y,c1,c2,x1,y1);

  x=n1;    y=0;      xy_rad(rc,dc,x,y,&alpha,&delta,1,0);
  toms2(alpha,c1,1); toms2(delta,c2,0); 
  rad_xy(rc,dc,alpha,delta,&x1,&y1,1,0);
  printf("%8.1lf %8.1lf-> %s %s-> %8.1lf %8.1lf\n",x,y,c1,c2,x1,y1);

  x=0;     y=0;      xy_rad(rc,dc,x,y,&alpha,&delta,1,0);
  toms2(alpha,c1,1); toms2(delta,c2,0); 
  rad_xy(rc,dc,alpha,delta,&x1,&y1,1,0);
  printf("%8.1lf %8.1lf-> %s %s-> %8.1lf %8.1lf\n",x,y,c1,c2,x1,y1);
/*
 i=0;     j=1;            pgbegin_(&i,"/xw",&j,&j,3);
 z0=17.;  z1=0.1;         pgpap_(&z0,&z1);
 z0=0.;   z1=n1;  z2=n2;  pgenv_(&z0,&z1,&z0,&z2,&i,&i);
*/
  printf("used file:\n");
  fp1=fopen("/line3/uband-data/cat/fits.pos","r");
l10:
  fgets(b,144,fp1); if(feof(fp1))goto l20;
  sscanf(&b[46],"%lf %lf %lf %lf %lf %lf %lf %lf",
  &aa[0],&dd[0],&aa[1],&dd[1],&aa[3],&dd[3],&aa[2],&dd[2]); 
//  if(b[40]!='4')goto l10;
//  for(i=0;i<4;i++)printf("%lf %lf ",aa[i],dd[i]); printf("\n");
  for(i=0;i<4;i++){
    rad_xy(rc,dc,aa[i],dd[i],&x1,&y1,1,0); // as Schmidt tel.
    xx[i]=x1; yy[i]=y1;
  } xx[4]=xx[0]; yy[4]=yy[0];
// i=5; pgline_(&i,xx,yy);
  for(i=0;i<4;i++){ ix[i]=xx[i]+0.5; iy[i]=yy[i]+0.5;  }
  i1=ix[0]; if(i1>ix[1])i1=ix[1];
  i2=ix[2]; if(i2<ix[3])i2=ix[3];
  j1=iy[0]; if(j1>iy[3])j1=iy[3];
  j2=iy[1]; if(j2<iy[2])j2=iy[2];
  if(j1<0)j1=0; if(j2>n2-1)j2=n2-1;
  if(i1<0)i1=0; if(i2>n1-1)i2=n1-1;
  if(j2-j1<10 || i2-i1<10)goto l10;
//     printf("%d %d    %d %d dj_di: %d %d\n",i1,j1,i2,j2,j2-j1,i2-i1);
  b[45]=0; shrinkv();          // b_filename. c[]_work, d[]_shrink_data
  i=ccdno-1; q1=q1d[i]; q2=4031.-q2u[i];
 printf("%s\n",b);
  for(j=j1;j<=j2;j++)for(i=i1;i<=i2;i++){
    x=i; y=j; // x+=0.5; y+=0.5;
    xy_rad(rc,dc,x,y,&alpha,&delta,1,0);  // big ix,iy
    rad_xy(s8[6],s8[7],alpha,delta,&x1,&y1,2,1);     // small x,y
    x=t8[0]*x1+t8[2]*y1+t8[4];
    y=t8[1]*x1+t8[3]*y1+t8[5];                       
    if(x<=0. || x>=4095)continue;         
    if(y<=q1 || y>=q2  )continue;          
    if(ccdno==4)if(x<=2047 || y>=2015)continue;   
    iy0=y;
    iy1=iy0+1;  
    wy=y-iy0;   zy=1.-wy;
    ix0=x;
    ix1=ix0+1;  
    wx=x-ix0;   zx=1.-wx;
    z0=zx*zy; z1=wx*zy; z2=zx*wy; z3=wx*wy;
    a[j*n1+i]+=c[iy0][ix0]*z0+c[iy0][ix1]*z1+c[iy1][ix0]*z2+c[iy1][ix1]*z3;
  }
  goto l10;
l20:
  fclose(fp1);
  puthead();
  fp=fopen(f1,"wb");
  fwrite(head,72,80,fp); swap4(a,n1*n2*4);
  fwrite(a,n1*n2,4,fp);  fclose(fp);
  printf("\nproduce file: %s    ok!\n",f1);
// pgend_();
}

shrinkv()
{
  int i,j,k,ii,jj,i1,j1;
  char h[72][80];
  fp=fopen(b,"rb");
  if(fp==0){ printf("%s not found!\n"); exit(0); }
  fread(h,72,80,fp);
  fread(c,4096*4032,4,fp); fclose(fp); swap4(c,4096*4032*4);
  k=indexpos(h,"CCD_NO: ",72);
  sscanf(&h[k][23],"%d",&ccdno);
  k=indexpos(h,"SKYADU  ",72);
  sscanf(&h[k][23],"%f",&white);
  k=indexpos(h,"A81     ",72);
  for(i=0;i<8;i++)sscanf(&h[k+i][10],"%le",&s8[i]);
  xytoad(s8,t8);  
} 

puthead()
{
  int i,k;
  double x;
  for(i=0;i<72;i++)for(k=0;k<80;k++)head[i][k]=32;
  strcpy(&head[ 0][0], "SIMPLE  =                    T ");
  strcpy(&head[ 1][0], "BITPIX  =                  -32 ");
  strcpy(&head[ 2][0], "NAXIS   =                    2 ");
  sprintf(&head[3], "NAXIS1  =                %5d ",n1);
  sprintf(&head[4], "NAXIS2  =                %5d ",n2);
  strcpy(&head[ 8][0], "DATE-OBS=     2010,10,02~12,09 ");
  strcpy(&head[ 9][0], "TIME-OBS=                      ");
  strcpy(&head[10][0], "EXPTIME =                  150 ");
  toms2(rc/cx,c1,1); toms2(dc/cy,c2,0);
  sprintf(&head[11],"RA      =        %s ", c1);
  sprintf(&head[12],"DEC     =         %s ",c2);
  strcpy(&head[13][0], "HA      =                      ");
  strcpy(&head[14][0], "EPOCH   =               2000.0 ");
  strcpy(&head[15][0], "FILTER  =                    u ");
  strcpy(&head[18][0], "INSTRUME=        BOK   90prime ");
  k=34;
  sprintf(&head[k++],"A81     = %20.12le ",a8[0]);
  sprintf(&head[k++],"A82     = %20.12le ",a8[1]);
  sprintf(&head[k++],"A83     = %20.12le ",a8[2]);
  sprintf(&head[k++],"A84     = %20.12le ",a8[3]);
  sprintf(&head[k++],"A85     = %20.12le ",a8[4]);
  sprintf(&head[k++],"A86     = %20.12le ",a8[5]);
  sprintf(&head[k++],"A87     = %20.12le ",a8[6]);
  sprintf(&head[k++],"A88     = %20.12le ",a8[7]);

  k=48;
  x=(a8[1]+a8[2])/(a8[0]-a8[3])/cy;       // BOKs
  sprintf(&head[k++],"CTYPE1  = 'RA---ARC'           ");
  sprintf(&head[k++],"CTYPE2  = 'DEC--ARC'           ");
  sprintf(&head[k++],"CROTA2  =  %16.8le    ",x);
  sprintf(&head[k++],"CRVAL1  =  %16.8le    ",rc/cy);
  sprintf(&head[k++],"CRVAL2  =  %16.8le    ",dc/cy);
  sprintf(&head[k++],"CDELT1  =  %16.8le    ",a8[0]/cy);
  sprintf(&head[k++],"CDELT2  =  %16.8le    ",a8[3]/cy);
//  sprintf(&head[k++],"CRPIX1  =  %16.8le    ",n1/2.);
//  sprintf(&head[k++],"CRPIX2  =  %16.8le    ",n2/2.);
  sprintf(&head[k++],"CRPIX1  =  %16.8le    ",(n1+1)/2.);
  sprintf(&head[k++],"CRPIX2  =  %16.8le    ",(n2+1)/2.);
// while nei_cha, use ix0,iy0, not ix1.iy1. record_down a[][], lost halt pixel 

  for(i=0;i<70;i++)if(head[i][0]!=32)head[i][31]='/';
  strcpy(&head[71][0], "END"); head[71][3]=32;
}

xytoad(x,a)
double *x,*a;
{
  double z;
  z=x[0]*x[3]-x[2]*x[1];
  a[0]= x[3]/z;
  a[1]=-x[1]/z;
  a[2]=-x[2]/z;
  a[3]= x[0]/z;
  a[4]= (x[2]*x[5]-x[4]*x[3])/z;
  a[5]= (x[4]*x[1]-x[0]*x[5])/z;
}

chs(a,y)
char *a; double *y;
{
  double x;
  char b[20];
  int i,j,k;
  strcpy(b,a); k=strlen(b);
  j=1; for(i=0;i<k;i++){ if(b[i]==32)continue; if(b[i]=='-')j=-1; break; }
  for(i=0;i<k;i++)if(b[i]<'.' || b[i]> '9')b[i]=32;
  i=k=0; x=0.; sscanf(b,"%d %d %lf",&i,&k,&x);
  x=k/60.+x/3600.+i;
  *y=x*j;
}

#define bok  48.
rad_xy(ac,dc,ra,de,x,y,type,k)
double ac,dc,ra,de,*x,*y;  int type,k;
{
  double xi,xn, tmp,rar,der;
  rar=ra*cx; der=de*cy;
/*                         // same discribe formula as following 3 lines
  double sd,cd,td,co;      // china bai_ke_quan_shu (astronmy) p.552
  sd=sin(dc);  cd=cos(dc);  td=tan(der);
  co=cos(rar-ac);  tmp=sd*td+cd*co;
  xi=sin(rar-ac)/tmp;
  xn=(cd*td-sd*co)/tmp;
*/
  tmp=atan(tan(der)/cos(rar-ac));
  xi=cos(tmp)*tan(rar-ac)/cos(tmp-dc);
  xn=tan(tmp-dc);
  if(type==1){     // schmidt TELESCOPE
    tmp=sqrt(xi*xi+xn*xn);
    if(tmp!=0.){
      tmp=atan(tmp)/tmp;
      xi*=tmp; xn*=tmp;
    }
  }
  if(type==2){    //  BOK telescope
    tmp=1.+bok*(xi*xi+xn*xn);
    xi*=tmp; xn*=tmp;
  }
  if(k==0){
    *x=b8[0]*xi+b8[2]*xn+b8[4];
    *y=b8[1]*xi+b8[3]*xn+b8[5];
  }
  if(k==1){ *x=xi; *y=xn;}
}

xy_rad(ac,dc,x,y,ra,de,type,k)
double ac,dc,x,y,*ra,*de;  int type,k;
{
  double xi,xn,tmp,tmp1;
  double x1,y1,r1,d1;
  int i;
  xi=a8[0]*x+a8[2]*y+a8[4];
  xn=a8[1]*x+a8[3]*y+a8[5];
  if(type==1){
    tmp=sqrt(xi*xi+xn*xn);
    if(tmp!=0.){
      tmp=tan(tmp)/tmp;
      xi*=tmp;  xn*=tmp;
    }
  }
  if(type==2){
    tmp=1.+bok*(xi*xi+xn*xn);

    for(i=0;i<4;i++){
      x1=xi/tmp; y1=xn/tmp;
      tmp=tan(dc);  tmp1=1.-y1*tmp;
      r1=atan(x1/cos(dc)/tmp1);
      d1=atan((y1+tmp)*cos(r1)/tmp1);
      r1+=ac;
      tmp=atan(tan(d1)/cos(r1-ac));
      x1=cos(tmp)*tan(r1-ac)/cos(tmp-dc);
      y1=tan(tmp-dc);
      tmp=1.+bok*(x1*x1+y1*y1);
//                                 printf("%lf\n",tmp);
    }
    xi/=tmp; xn/=tmp;
  }
  tmp=tan(dc);  tmp1=1.-xn*tmp;
  *ra=atan(xi/cos(dc)/tmp1);
  *de=atan((xn+tmp)*cos(*ra)/tmp1);
  *ra+=ac;
  tmp=(*de)-dc; if(tmp<0.)tmp=-tmp;
  if(tmp>2.){  *de=-(*de);  *ra+=pi; }
  if(*ra<0.)*ra+=(pi+pi);
  if(k==0){ *ra/=cx; *de/=cy; }
}

toms2(aa,c,k)
double aa;           // input
char c[14];          // output
int k;               // input
/*
 hour or degree to char_line
 if k=1 (hour case) in **:**:**.***
    k=0 (degree       -**:**:**.**
*/
{
  int  i,j;
  double a,b,x;
  a=aa;
  if(k==1){ if(a>24.)a-=24.; if(a<0.)a+=24.; }
  c[0]=32; if(a<0.){ a=-a; c[0]='-'; }
  i=a; b=(a-i)*60.;
  j=b; x=(b-j)*60.;
  if(j==60){ j=0; i++; }
  if(k==1)sprintf(&c[1],"%2.2d:%2.2d:%6.3f",i,j,x);
  if(k==0)sprintf(&c[1],"%2.2d:%2.2d:%5.2f",i,j,x);
  if(c[7]==32)c[7]='0';   if(c[8]==32)c[8]='0';
  if(c[7]=='6'){
    c[7]='0'; j++; sprintf(&c[4],"%2.2d",j);
    if(c[4]=='6'){
      c[4]='0'; i++; sprintf(&c[1],"%2.2d",i);
    }
    c[3]=':'; c[6]=':';
  }
}