#include <stdio.h>
#include <math.h>
#define n1   28400            // X,n1
#define n2    5000            // Y,n2
#define scale 15.0            // extract from *.fit in step of 10 aecsec/pixel
#define f1  "00:45:00"        // R.A. center
#define f2  "-01:00:0"        // DEC. center
#define v       32            // shrink *.fits to 16*0.45=14.5 arcsec/pixel
#define m1      4096/v
#define m2      4032/v
char   head[72][80],b[146];
float  a[n2][n1];             // a[n2][n1]  a[y][x]  
float  c[4032][4096],d[m2][m1];
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],ccdno,idate,ip=0,pp;
FILE   *fp,*fp1;
int    q1d[2][4],q2u[2][4],q3l[2][4],q4r[2][4];
unsigned char  bb[n2*n1];
main(int ac, char **av)
{
  int    i,j,k,i1,i2,j1,j2;
  double x,y,alpha,delta,x1,y1;
  int    ix0,ix1,iy0,iy1;
  float  z0,z1,z2,z3,zz;
  float  zx,zy,wx,wy,flux;
  float  z2016_1,z2048,z4032_1,z4096_1;
  int    p1,p2,q1,q2;
  printf("\n\t v_72 [!]\n\n");
  if(ac<2)pp=1; else pp=0;
  fp=fopen("/vega2/rhbin/ubandlines.dat","r");
  if(fp==0){ printf("ubandlines.dat not found!\n"); exit(0); }
  for(i=0;i<4;i++){ 
l00:   fgets(b,80,fp);  if(b[0]=='#')goto l00;
       sscanf(b,"%d %d %d %d",&q1d[0][i],&q2u[0][i],&q3l[0][i],&q4r[0][i]);}
  for(i=0;i<4;i++){ 
l01:   fgets(b,80,fp);  if(b[0]=='#')goto l01;
       sscanf(b,"%d %d %d %d",&q1d[1][i],&q2u[1][i],&q3l[1][i],&q4r[1][i]);}
  fclose(fp);

  z2016_1=2016./v-1.;
  z2048  =2048./v;
  z4032_1=4032./v-1.;
  z4096_1=4096./v-1.;
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  for(j=0;j<n2;j++)for(i=0;i<n1;i++)a[j][i]=0.;
  printf("center:  %s %s\n",f1,f2);
  chs(f1,&rc); chs(f2,&dc); 
  rc*=cx; dc*=cy; a8[6]=rc; a8[7]=dc;  
  hs2=(3600.*180./pi)/scale;                 
  b8[0]=-hs2; b8[1]=b8[2]=0; b8[3]=-b8[0]; b8[4]=n1/2; b8[5]=n2/2;
  xytoad(b8,a8);
   
  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);
   
  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=1;     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=1;      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=1;     y=1;      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;              if(pp)pgbegin_(&i,"/xw",&j,&j,3);
 z0=12.;  z1=.25;           if(pp)pgpap_(&z0,&z1);
 zx=0; z1=n1; zy=0; z2=n2;  if(pp)pgenv_(&zx,&z1,&zy,&z2,&i,&i);

  fp1=fopen("/line3/uband-data/cat/fits.pos9","r");
l10:
  fgets(b,146,fp1); if(feof(fp1))goto l20;
  if(b[0]=='#')goto l10;
  ip++; 
//      if(ip%10)goto l10;

// if(ip<3670)goto l10;
// if(ip>3750)goto l20;

  flux=1;
  sscanf(&b[46],"%lf %lf %lf %lf %lf %lf %lf %lf %f",
  &aa[0],&dd[0],&aa[1],&dd[1],&aa[3],&dd[3],&aa[2],&dd[2],&flux); 
  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; if(pp)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;
  b[45]=0; shrinkv();          // b_filename. c[]_work, d[]_shrink_data
  i=ccdno-1; 
  sscanf(&b[30],"%d",&k);
  idate=0; if(k>5700)idate=1;
  printf("%d: %s %d%7.3f%7.3f%7.3f\n",ip,&b[29],idate,aa[0],dd[0],flux);
  p1=(q3l[idate][i]+31)/v; p2=(q4r[idate][i]+31)/v; p2=z4096_1-p2;
  q1=(q1d[idate][i]+31)/v; q2=(q2u[idate][i]+31)/v; q2=z4032_1-q2;
  for(j=j1;j<=j2;j++)for(i=i1;i<=i2;i++){
    x=i; y=j;              
    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];                       // use one side data
    x/=v; y/=v;
    if(x<=p1 || x>=p2 )continue;
    if(y<=q1 || y>=q2 )continue;
    if(idate==0 && ccdno==4)if(x<=z2048 || y>=z2016_1)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][i]+=d[iy0][ix0]*z0+d[iy0][ix1]*z1+d[iy1][ix0]*z2+d[iy1][ix1]*z3;
  }
  goto l10;
l20:
  fclose(fp1);
  puthead();
  head[1][27]=32; head[1][28]=' '; head[1][29]='8';
  i1=0;
  for(j=0;j<n2;j++)for(i=0;i<n1;i++){ 
    k=a[j][i]+0.5; if(k>255)k=255; bb[i1++]=k;
  }
  fp=fopen("uband_10370.fit","wb"); 
  fwrite(head,72,80,fp); fwrite(bb,n1,n2,fp); fclose(fp);
  printf("\nproduce file: uband_10370.fit    ok!\n");
 if(pp)pgend_();
}

shrinkv()
{
  int i,j,k,ii,jj,i1,j1;
  float x,v2,white,ww;
  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);  white+=10.;
  k=indexpos(h,"A81     ",72);
  for(i=0;i<8;i++)sscanf(&h[k+i][10],"%le",&s8[i]);
  xytoad(s8,t8);  v2=v*v; k++;
  ww=white*1.01;
  for(j=0;j<m2;j++)for(i=0;i<m1;i++){
    ii=i*v; jj=j*v; x=0.;
    for(j1=jj;j1<jj+v;j1++)for(i1=ii;i1<ii+v;i1++)x+=c[j1][i1];
    x/=v2; x-=ww; x/=2.; if(x<0.)x=0.; d[j][i]=x+1.;
  }
} 

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=        2011,09,18~24 ");
  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+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]=':';
  }
}