#include <stdio.h>
#include <stdlib.h>
#include <string.h>


void drotg_(double *DA, double *DB, double *C, double *S);
void drot_(const int *N, double*DX, const int *INCX, double *DY, 
	const int *INCY, const double *C, const double *S);

#define index2D(i,j,LD1) i + (j*LD1)
#define index3D(i,j,k,LD1,LD2) i + (j*LD1) + (k*LD1*LD2)


int MINVAL(int s, int*a)
{
   int v, e;

   v = a[0];

   for ( e = 0; e < s; e++ )
   {
       if ( v > a[e] ) v = a[e];
   
   }
   
   return(v);
}

int MAXVAL(int s, int*a)
{
   int v, e;

   v = a[0];

   for ( e = 0; e < s; e++ )
   {
       if ( v < a[e] ) v = a[e];
   
   }
   
   return(v);
}
 
void IntData2pgm(int s1, int s2, int* idata, char* name) 
{
   /*
      Simple subroutine to dump integer data in a PGM format
   */
   int i, j;
   FILE* ouni;
   int  vp, vs;
   float  rmin, rmax;
   char  fname[80];
   /*
      Write on unit 700 with  PGM format
   */
   strcpy(fname,name);
   strcat(fname,".pgm\0");
   ouni = fopen(fname,"w");
   if ( ! ouni ) {
      fprintf(stderr,"!!!! Error write access to file %s\n",fname);
   }
   /*  Magic code */
   fprintf(ouni,"P2\n");
   /*  Dimensions */
   fprintf(ouni,"%d %d\n",s1, s2);
   /*  Maximum value */
   fprintf(ouni,"255\n");
   /*  Values from 0 to 255 */
   rmin = MINVAL(s1*s2,idata); rmax = MAXVAL(s1*s2,idata);
   vs = 0;
   for ( i = 0; i < s1; i++ ) {
      for ( j = 0; j < s2; j++ ) {
         vp = (int) ( (idata[i+(j*s1)] - rmin) * 255.0 / (rmax - rmin) );
         vs++;
         fprintf(ouni,"%4d ", vp);
         if ( vs >= 10 ) {
            fprintf(ouni," \n");
            vs = 0;
         }
      }
      fprintf(ouni," ");
      vs = 0;
   }
   fclose(ouni);
   return;
}

   
void DrawFigure(int ld, int n, int *plane, 
	     int lds, int ns, double *square, double opaq)
{
	/*
	Given a matrix plane[ld,n] representing a rectangle of points in the 2D plane
	points are set with opaq opacity coefficient according to definitions
	in square[lds,ns,2] matrix
	*/
	int i, j, ix, iy, pv;
	
	pv = (int)( (double)256*opaq )-1; 
	if ( pv > 255 ) pv = 255;
   for ( i = 0; i < ns; i++ ) {
      for ( j = 0; j < ns; j++ ) {
         if ( i == 0 || i == ns-1 || j == 0 || j == ns-1 ) {
            ix = square[index3D(i,j,0,lds,ns)]+n/2; 
            iy = square[index3D(i,j,1,lds,ns)]+n/2;
            if ( ix > 0 && ix < n ) {
            	if ( iy > 0 && iy < n ) {
            		plane[index2D(ix,iy,ld)] = pv;	// Draw figure into plane
            	}
            }
         }
      }
   }
   return;
}

int main(  int argc, char *argv[] )
{ 

//
//   Use BLAS subroutines DROTG and DROT to generate a rotation and calculate the 
//      rotation of a squared figure 
//
   int *plane;   // Drawing 2D rectangle
   double  *square;   // Squared figure
   int ld, n, lds, ns, rc, i, j, ix, iy, ii, jj, dims, one;
   double da, db, c, s, opaq;
   
   n = 0;
   printf("Plane dimension?\n");
   fscanf(stdin,"%d",&n);
   if ( n < 200 || n > 1000 ) n = 1000;
   printf("dimension = %d\n",n);
   ld = n;
   plane = (int*) malloc(sizeof((int)1)*ld*n);
   if ( plane == NULL ) return(-1);

   ns = n/4; lds = ns;
   square = (double*) malloc(sizeof((double)1.0)*lds*ns*2);
   if ( square == NULL ) return(-1);
   
   for ( i = 0; i < ld*n; i++ ) plane[i] = 0.0; 	// Initialize plane
   //
   // Define squared figure ...
      for ( i = 0; i < ns; i++ ) {
      	  for ( j = 0; j < ns; j++ ) {
      	  	  if ( i == 0 || i == ns-1 || j == 0 || j == ns-1 ) {
      	  	  	  square[index3D(i,j,0,lds,ns)] = -ns/2 + i;
      	  	  	  square[index3D(i,j,1,lds,ns)] = -ns/2 + j;
      	  	  }
      	  }
      }

   // ... and trace it in plane
   opaq = 0.5;
   DrawFigure(ld,n,plane,lds,ns,square,opaq);
   
   printf("Define rotation: da, db ?\n");
   fscanf(stdin,"%lf, %lf",&da,&db);
//
   drotg_(&da,&db,&c,&s);	// get parameters C, S
   printf("Rotation parameters: %lf, %lf\n",c,s);
   dims = ns*ns;
   one = 1;
   drot_(&dims,&square[index3D(0,0,0,lds,ns)],&one,
   	       &square[index3D(0,0,1,lds,ns)],&one,&c,&s);
//       
//	Draw new figure in plane
   opaq = 1.0;
   DrawFigure(ld,n,plane,lds,ns,square,opaq);
//
   IntData2pgm(n,n,plane,"plane");
//   
   return(0);
}