# include "data.h"

int rotate (int argc, char **argv) {
  
  int i, j, NX, NY, X, Y, Lx, Ly, N, newCenter;
  float *in_buff, *out_buff, *c;
  double angle, CosAngle, SinAngle, Xo, Yo, dX, dY, fx, fy, x, y, X1, Y1;
  double pc11, pc12, pc21, pc22, PC11, PC12, PC21, PC22;
  Buffer *buf;

  Xo = 0;
  Yo = 0;
  if (N = get_argument (argc, argv, "-center")) {
    newCenter = TRUE;
    remove_argument (N, &argc, argv);
    Xo  = atof(argv[N]);
    remove_argument (N, &argc, argv);
    Yo  = atof(argv[N]);
    remove_argument (N, &argc, argv);
  }
  else 
    newCenter  = FALSE;

  if (argc != 3) {
    fprintf (stderr, "USAGE: rotate <buffer> <angle>\n");
    return (FALSE);
  }
  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);

  if ((atof (argv[1]) < -180) || (atof (argv[1]) > 180)) {
    fprintf (stderr, "valid rotate angle between -180 and +180 degrees\n");
    return (FALSE);
  }

  /* save starting values */
  NX = buf[0].header.Naxis[0];
  NY = buf[0].header.Naxis[1];
  in_buff = (float *) buf[0].matrix.buffer;  /* don't lose reference */

  if (!strcasecmp (argv[2], "LEFT") || (atof (argv[2]) == -90)) {
    buf[0].header.Naxis[0] = NY;
    buf[0].header.Naxis[1] = NX;
    fits_modify (&buf[0].header, "NAXIS1", "%d", 1, NY);
    fits_modify (&buf[0].header, "NAXIS2", "%d", 1, NX);
    fits_print  (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!");
    fits_create_matrix (&buf[0].header, &buf[0].matrix);
    out_buff = (float *)buf[0].matrix.buffer;
    for (i = NX - 1; i > -1; i--) {
      for (j = 0; j < NY; j++, out_buff++) {
	*out_buff = in_buff[i + j*NX];
      }
    }
    /* fix reference pixel */
    fits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo);
    fits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo);
    X1 = Yo;
    Y1 = NX - Xo;
    fits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1);
    fits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1);
    
    /* fix rotate matrix */
    fits_scan (&buf[0].header, "PC001001", "%lf", 1, &pc11);
    fits_scan (&buf[0].header, "PC001002", "%lf", 1, &pc12);
    fits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21);
    fits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22);
    PC11 = pc21;
    PC12 = pc22;
    PC21 = -pc11;
    PC22 = -pc12;
    fits_modify (&buf[0].header, "PC001001", "%le", 1, PC11);
    fits_modify (&buf[0].header, "PC001002", "%le", 1, PC12);
    fits_modify (&buf[0].header, "PC002001", "%le", 1, PC21);
    fits_modify (&buf[0].header, "PC002002", "%le", 1, PC22);

    free (in_buff);
    return (TRUE);
  }

  if (!strcasecmp (argv[2], "RIGHT") || (atof (argv[2]) == 90)) {
    buf[0].header.Naxis[0] = NY;
    buf[0].header.Naxis[1] = NX;
    fits_modify (&buf[0].header, "NAXIS1", "%d", 1, NY);
    fits_modify (&buf[0].header, "NAXIS2", "%d", 1, NX);
    fits_print  (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!");
    fits_create_matrix (&buf[0].header, &buf[0].matrix);
    out_buff = (float *)buf[0].matrix.buffer;
    for (i = 0; i < NX; i++) {
      for (j = NY - 1; j > -1; j--, out_buff++) {
	*out_buff = in_buff[i + j*NX];
      }
    }
    /* fix reference pixel */
    fits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo);
    fits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo);
    X1 = NY - Yo;
    Y1 = Xo;
    fprintf (stderr, "%f %f -> %f %f\n", Xo, Yo, X1, Y1);
    fits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1);
    fits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1);
    
    /* fix rotate matrix */
    fits_scan (&buf[0].header, "PC001001", "%lf", 1, &pc11);
    fits_scan (&buf[0].header, "PC001002", "%lf", 1, &pc12);
    fits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21);
    fits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22);
    PC11 = -pc21;
    PC12 = -pc22;
    PC21 = pc11;
    PC22 = pc12;
    fprintf (stderr, "%f %f  ->  %f %f\n", pc11, pc12, PC11, PC12);
    fprintf (stderr, "%f %f  ->  %f %f\n", pc21, pc22, PC21, PC22);
    fits_modify (&buf[0].header, "PC001001", "%le", 1, PC11);
    fits_modify (&buf[0].header, "PC001002", "%le", 1, PC12);
    fits_modify (&buf[0].header, "PC002001", "%le", 1, PC21);
    fits_modify (&buf[0].header, "PC002002", "%le", 1, PC22);

    free (in_buff);
    return (TRUE);
  }

  if (!strcasecmp (argv[2], "UPSIDE") || (atof (argv[2]) == -180) || (atof (argv[2]) == 180)) {
    fits_print  (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!");
    fits_create_matrix (&buf[0].header, &buf[0].matrix);
    out_buff = (float *)buf[0].matrix.buffer;
    for (j = NY - 1; j > -1; j--) {
      for (i = NX - 1; i > -1; i--, out_buff++) {
	*out_buff = in_buff[i + j*NX];
      }
    }
    /* fix reference pixel */
    fits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo);
    fits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo);
    X1 = NX - Xo;
    Y1 = NY - Yo;
    fits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1);
    fits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1);
    
    /* fix rotate matrix */
    fits_scan (&buf[0].header, "PC001001", "%lf", 1, &pc11);
    fits_scan (&buf[0].header, "PC001002", "%lf", 1, &pc12);
    fits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21);
    fits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22);
    PC11 = -pc11;
    PC12 = -pc12;
    PC21 = -pc21;
    PC22 = -pc22;
    fits_modify (&buf[0].header, "PC001001", "%le", 1, PC11);
    fits_modify (&buf[0].header, "PC001002", "%le", 1, PC12);
    fits_modify (&buf[0].header, "PC002001", "%le", 1, PC21);
    fits_modify (&buf[0].header, "PC002002", "%le", 1, PC22);

    free (in_buff);
    return (TRUE);
  }

  if (!strcasecmp (argv[2], "FLIPY")) {
    fits_print  (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!");
    fits_create_matrix (&buf[0].header, &buf[0].matrix);
    out_buff = (float *)buf[0].matrix.buffer;
    for (j = NY - 1; j > -1; j--) {
      for (i = 0; i < NX; i++, out_buff++) {
	*out_buff = in_buff[i + j*NX];
      }
    }
    /* fix reference pixel */
    fits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo);
    fits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo);
    X1 = Xo;
    Y1 = NY - Yo;
    fits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1);
    fits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1);
    
    /* fix rotate matrix */
    fits_scan (&buf[0].header, "PC001001", "%lf", 1, &pc11);
    fits_scan (&buf[0].header, "PC001002", "%lf", 1, &pc12);
    fits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21);
    fits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22);
    PC11 = pc11;
    PC12 = -pc12;
    PC21 = pc21;
    PC22 = -pc22;
    fits_modify (&buf[0].header, "PC001001", "%le", 1, PC11);
    fits_modify (&buf[0].header, "PC001002", "%le", 1, PC12);
    fits_modify (&buf[0].header, "PC002001", "%le", 1, PC21);
    fits_modify (&buf[0].header, "PC002002", "%le", 1, PC22);

    free (in_buff);
    return (TRUE);
  }

  if (!strcasecmp (argv[2], "FLIPX")) {
    fits_print  (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!");
    fits_create_matrix (&buf[0].header, &buf[0].matrix);
    out_buff = (float *)buf[0].matrix.buffer;
    for (j = 0; j < NY; j++) {
      for (i = NX - 1; i > -1; i--, out_buff++) {
	*out_buff = in_buff[i + j*NX];
      }
    }
    /* fix reference pixel */
    fits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo);
    fits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo);
    X1 = NX - Xo;
    Y1 = Yo;
    fits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1);
    fits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1);
    
    /* fix rotate matrix */
    fits_scan (&buf[0].header, "PC001001", "%lf", 1, &pc11);
    fits_scan (&buf[0].header, "PC001002", "%lf", 1, &pc12);
    fits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21);
    fits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22);
    PC11 = -pc11;
    PC12 = pc12;
    PC21 = -pc21;
    PC22 = pc22;
    fits_modify (&buf[0].header, "PC001001", "%le", 1, PC11);
    fits_modify (&buf[0].header, "PC001002", "%le", 1, PC12);
    fits_modify (&buf[0].header, "PC002001", "%le", 1, PC21);
    fits_modify (&buf[0].header, "PC002002", "%le", 1, PC22);

    free (in_buff);
    return (TRUE);
  }

  angle = atof (argv[2]);
  CosAngle = cos (angle*RAD_DEG);
  SinAngle = sin (angle*RAD_DEG);
  
  fprintf (stderr, "rotating: %f %f %f\n", angle, CosAngle, SinAngle);

  Lx = NX*fabs(CosAngle) + NY*fabs(SinAngle);
  Ly = NX*fabs(SinAngle) + NY*fabs(CosAngle);
  dX = MAX(0,NY*SinAngle);
  dY = MAX(0,-NX*SinAngle);
  /*
  fprintf (stderr, "%f %f  -->  ", Xo, Yo);
  X1 = Xo*CosAngle - Yo*SinAngle + dX;
  Y1 =  Xo*SinAngle + Yo*CosAngle + dY;
  fprintf (stderr, "%f %f\n", X1, Y1);
  */

  /* fix reference pixel */
  fits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo);
  fits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo);
  /*
  X1 = (Xo - dX)*CosAngle - (Yo - dY)*SinAngle;
  Y1 = (dX - Xo)*SinAngle + (Yo - dY)*CosAngle;
  */
  X1 = Xo*CosAngle - Yo*SinAngle + dX;
  Y1 =  Xo*SinAngle + Yo*CosAngle + dY;
  fits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1);
  fits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1);

  /* fix rotate matrix */
  fits_scan (&buf[0].header, "PC001001", "%lf", 1, &pc11);
  fits_scan (&buf[0].header, "PC001002", "%lf", 1, &pc12);
  fits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21);
  fits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22);
  PC11 = pc11*CosAngle - pc21*SinAngle;
  PC12 = pc12*CosAngle - pc22*SinAngle;
  PC21 = pc21*CosAngle + pc11*SinAngle;
  PC22 = pc22*CosAngle + pc12*SinAngle;
  fits_modify (&buf[0].header, "PC001001", "%le", 1, PC11);
  fits_modify (&buf[0].header, "PC001002", "%le", 1, PC12);
  fits_modify (&buf[0].header, "PC002001", "%le", 1, PC21);
  fits_modify (&buf[0].header, "PC002002", "%le", 1, PC22);

  buf[0].header.Naxis[0] = Lx;
  buf[0].header.Naxis[1] = Ly;
  fits_modify (&buf[0].header, "NAXIS1", "%d", 1, Lx);
  fits_modify (&buf[0].header, "NAXIS2", "%d", 1, Ly);
  fits_create_matrix (&buf[0].header, &buf[0].matrix);
  fits_print  (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!");
  out_buff = (float *)buf[0].matrix.buffer;
  for (j = 0; j < Ly; j++) {
    for (i = 0; i < Lx; i++, out_buff++) {

      x = (i - dX)*CosAngle + (j - dY)*SinAngle;
      y = (dX - i)*SinAngle + (j - dY)*CosAngle;
      X = (int) x;
      Y = (int) y;

      if (X < 0) continue;
      if (X >= NX - 1) continue;
      if (Y < 0) continue;
      if (Y >= NY - 1) continue;

      c = &in_buff[X * NX*Y];
      fx = x - X;
      fy = y - Y;
      *out_buff = (c[0]*(1-fx) + c[1]*fx)*(1-fy) + (c[NX+1]*fx + c[NX]*(1-fx))*fy;
    }
  }
  free (in_buff);
  return (TRUE);

}
