/* Rotate and display a radiation plot correctly from any angle */ 
/* this uses heap-sort, so array begins at 1 and ends at N */
/* Uses all three Euler angles. for physics 302  */
/* gcc -o radiation radiation.c -L/usr/X11R6/lib -lplot -lXaw -lXmu -lXt -lSM  -lICE -lXext -lX11 -lm */
/* Jeffrey R. Schmidt 2000 */
/* Read through all of the comments if you want to figure out */
/* how this program works and make useful modifications */
/* capable of X, PS or GIF output specified on command line */
/* can cut solid figure in half two ways */

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<plot.h>
#define N 2048    /* this is number of plaquettes to be read */
#define Np 32 /* number of azimuthal divisions */
#define Nt 65 /* number of latitude divisions */
#define PI 3.1415926

struct plaquette {
    /* to learn about this, read my 499 handout on graphics */
  float vertex[4][3];  /* x and y coords of basepoint and other 3 points */
  float center[3];     /* x and y coordinates of its center */
  float norml[3];     /* normal vector to plaquette face */
  float face;
  float vdist;         /* distance from vantage point */
  int clr;             /* grayscale shade             */
} pl[N+2],tempPL;


float dist(float x1, float y1, float z1, float x2, float y2, float z2);
void getplaquetteattribs( struct plaquette *Plaq, float u, float v, float w);
float function(float x, float y);
int mod(int A, int B);

float dphi,dtheta,Stheta,Sphi,Srad, vantx,vanty,vantz,R;
float a11,a12,a13,a21,a22,a23,a31,a32,a33;
float Z[Np*Nt+2],X[Np*Nt+2],Y[Np*Nt+2];
int n,m,CL,i,j,l,ir,k,u,u1,u2,u3;
float dimensionx,vantagex, vantagey,vantagez;/* position of a vantage point */ 
float theta,phi,psi,tmpx,tmpy,tmpz,xoffset,yoffset,resize;
int pl_handle,halfplot,displayt;
FILE *fptr;

main(int argc, char *argv[]){
dtheta=PI/(float)Nt;
dphi=2.0*PI/(float)Np;
vantagex=0.0;
vantagey=0.0;
vantagez=5.0;
 if(argc < 9 || argc > 9){
     printf("./radiation theta phi psi x_0 y_0 resize cut_param display\n");
     printf("Three Euler angles (radians) for view perspective.\n");
     printf("x_0,y_0  move the origin left/right, up/down\n");
     printf("Resize factor zoomes in or out, = 1.0 default\n");
     printf("Cut_param=0 for full plot, 1 for phi-cut, 2 for theta-cut\n");
     printf("Display=0 for X, 1 for PS, 2 for GIF\n");
     printf("THERE IS NO ERROR CHECKING; just a crash if your parameters suck\n");
     exit(1);
 }
 /* Euler angles to rotate viewpoint */
theta=atof(argv[1]);
phi=atof(argv[2]);
psi=atof(argv[3]);
/* move the origin left or right */
xoffset=atof(argv[4]);
yoffset=atof(argv[5]);
resize=atof(argv[6]);
halfplot=atoi(argv[7]);
displayt=atoi(argv[8]);


 if(halfplot == 1 ) dphi=dphi/2.0;
 if(halfplot == 2 ) dtheta=dtheta/2.0;
 

  /* create a spherical grid and get the scale */ 
  X[0]=0.0;Y[0]=0.0;Z[0]=function(0.0,0.0); /* north pole */
  dimensionx=0.3; /* default size of affine space */
  for(n=1;n<=Nt;n++){
    for(m=0;m<=Np;m++){
      Sphi=(float)m*dphi;
      Stheta=(float)n*dtheta;
      /* now remap the surface */
      R=function(Sphi,Stheta);
      if(R>dimensionx) dimensionx=R; /* this is the affine-space resize */
      X[Np*(n-1)+m+1]=R*sin((float)n*dtheta)*cos((float)m*dphi);
      Y[Np*(n-1)+m+1]=R*sin((float)n*dtheta)*sin((float)m*dphi);
      Z[Np*(n-1)+m+1]=R*cos((float)n*dtheta);}}
  /* create rotation matrix to transform to new vantage point */
 a11=cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi);
 a12=cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi);
 a13=sin(psi)*sin(theta);
 a21=-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi);
 a22=-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi);
 a23=cos(psi)*sin(theta);
 a31=sin(theta)*sin(phi);
 a32=-sin(theta)*cos(phi);
 a33=cos(theta);

for(n=1;n<Nt;n++){
for(m=0;m<Np;m++){
      /* this is where plaquette vertices are assigned */
      if(halfplot==0){
u=Np*(n-1)+m+1;
u1=Np*(n-1)+mod(m+1,Np)+1;
u2=Np*n+mod(m+1,Np)+1;
u3=Np*n+m+1;}
      if (halfplot==1){
	  /* one splitting */
u=Np*(n-1)+m+1;
u1=Np*(n-1)+mod(m+1,Np)+1;
u2=Np*n+mod(m+1,Np)+1;
u3=Np*n+m+1;
 if(m==Np-1 || m==0){
     u1=u;
     u2=u;
     u3=u;}
      }
      if(halfplot==2){
	  /* the other type of splitting */
u=Np*(n-1)+m+1;
u1=Np*(n-1)+mod(m+1,Np)+1;
u2=Np*n+mod(m+1,Np)+1;
u3=Np*n+m+1;
 if(n==Nt-1 || n==1){
     u1=u;
     u2=u;
     u3=u;}
      }

 pl[u].vertex[0][0]=X[u];
 pl[u].vertex[0][1]=Y[u];
 pl[u].vertex[0][2]=Z[u];
 pl[u].vertex[1][0]=X[u1];
 pl[u].vertex[1][1]=Y[u1];
 pl[u].vertex[1][2]=Z[u1];
 pl[u].vertex[2][0]=X[u2];
 pl[u].vertex[2][1]=Y[u2];
 pl[u].vertex[2][2]=Z[u2];
 pl[u].vertex[3][0]=X[u3];
 pl[u].vertex[3][1]=Y[u3];
 pl[u].vertex[3][2]=Z[u3];
  }}
/* the affine space is now about twice as wide and high as the largest */
/* dimension of the plot, so cut it down, remember you can translate */
/* the origin to re-center the figure */
 dimensionx=resize*0.8*dimensionx;
 printf("Plot will be %f X %f\n", 2.0*dimensionx, 2.0*dimensionx);

/* rotate the plaquette points */
 for(n=1;n<=N;n++){
   for(m=0;m<4;m++){
tmpx=a11*pl[n].vertex[m][0]+a12*pl[n].vertex[m][1]+a13*pl[n].vertex[m][2];
tmpy=a21*pl[n].vertex[m][0]+a22*pl[n].vertex[m][1]+a23*pl[n].vertex[m][2];
tmpz=a31*pl[n].vertex[m][0]+a32*pl[n].vertex[m][1]+a33*pl[n].vertex[m][2];
pl[n].vertex[m][0]=tmpx;
pl[n].vertex[m][1]=tmpy;
pl[n].vertex[m][2]=tmpz;
}}

 /* move the plot right or left */
 for(n=1;n<=N;n++){
   for(m=0;m<4;m++){
pl[n].vertex[m][0]=pl[n].vertex[m][0]-xoffset;
pl[n].vertex[m][1]=pl[n].vertex[m][1]-yoffset;
}}

 /* now compute positions of centers and normals */
 for(n=1;n<=N;n++){
pl[n].center[0]=0.25*(pl[n].vertex[0][0]+pl[n].vertex[1][0]+ pl[n].vertex[2][0]+pl[n].vertex[3][0]);
pl[n].center[1]=0.25*(pl[n].vertex[0][1]+pl[n].vertex[1][1]+ pl[n].vertex[2][1]+pl[n].vertex[3][1]);
pl[n].center[2]=0.25*(pl[n].vertex[0][2]+pl[n].vertex[1][2]+ pl[n].vertex[2][2]+pl[n].vertex[3][2]);
 pl[n].vdist=dist(vantagex,vantagey,vantagez,pl[n].center[0],pl[n].center[1],pl[n].center[2]);
getplaquetteattribs(&pl[n],1.0,0.3,0.2);


 }

/* now heap sort according to distance from vantage point , least comes first */

	l=N/2+1;
	ir=N;
	do {
		if (l > 1) {
                        l=l-1;
			tempPL=pl[l];
		} else {
			tempPL=pl[ir];
			pl[ir]=pl[1];
                        ir=ir-1;
                        if(ir==1){
				pl[1]=tempPL;
				break;
			}
		}
		i=l;
		j=l+l;
		do {
			if (j<ir && pl[j].vdist < pl[j+1].vdist) j=j+1;
			if (tempPL.vdist < pl[j].vdist) {
				pl[i]=pl[j];
				i=j;
				j=2*j;
			} else j=ir+1;
		}while (j <= ir);
		pl[i]=tempPL;

	}while(l>=1);
	if(displayt==0){
	    pl_handle = pl_newpl ("X", stdin, stdout, stderr);}
	if(displayt==1){
	    fptr=fopen("out.ps","w+");
	    pl_handle = pl_newpl ("PS", stdin, fptr, stderr);}
	if(displayt==2){
  fptr=fopen("out.gif","w+");
 pl_parampl ("BITMAPSIZE", "350x350");
 pl_parampl ("BG_COLOR", "white");
 pl_parampl ("TRANSPARENT_COLOR", "white");
 pl_handle = pl_newpl ("GIF", stdin, fptr, stderr);}
  
 

 
  pl_selectpl (pl_handle);
  if (pl_openpl () < 0)         /* open Plotter */
    {
      fprintf (stderr, "Couldn't open Plotter\n");
      return 1;
    }
pl_fspace (-dimensionx,-dimensionx,dimensionx, dimensionx);
pl_flinewidth (0.0005);
pl_erase();
pl_filltype(1);
 for(j=0;j<N;j++){
n=N-j;             /* reverse heap-sort ordering */

/* this is where you can shade the plaquettes according to lighting vantage */
pl_color(pl[n].clr,pl[n].clr,pl[n].clr); 
pl_fillcolor(pl[n].clr,pl[n].clr,pl[n].clr); 



pl_fmove(pl[n].vertex[0][0],pl[n].vertex[0][1]);
pl_fline(pl[n].vertex[0][0],pl[n].vertex[0][1],pl[n].vertex[1][0],pl[n].vertex[1][1]);
pl_fline(pl[n].vertex[1][0],pl[n].vertex[1][1],pl[n].vertex[2][0],pl[n].vertex[2][1]);
pl_fline(pl[n].vertex[2][0],pl[n].vertex[2][1],pl[n].vertex[3][0],pl[n].vertex[3][1]);
 pl_fline(pl[n].vertex[3][0],pl[n].vertex[3][1],pl[n].vertex[0][0],pl[n].vertex[0][1]);
 pl_endpath(); /* this is an important step ! */
 pl_colorname("black");
 pl_fline(pl[n].vertex[2][0],pl[n].vertex[2][1],pl[n].vertex[3][0],pl[n].vertex[3][1]);
pl_fline(pl[n].vertex[0][0],pl[n].vertex[0][1],pl[n].vertex[1][0],pl[n].vertex[1][1]);
 pl_fline(pl[n].vertex[3][0],pl[n].vertex[3][1],pl[n].vertex[0][0],pl[n].vertex[0][1]);
pl_fline(pl[n].vertex[1][0],pl[n].vertex[1][1],pl[n].vertex[2][0],pl[n].vertex[2][1]);
}

 /* delete the plotters, and close any open file streams*/

  pl_closepl ();
  pl_selectpl (0);
  pl_deletepl (pl_handle);
 if(displayt==1 || displayt==2) fclose(fptr);
}

/* requisite subroutines */
float dist(float x1, float y1, float z1, float x2, float y2, float z2){
  /* computes squared point to point distance */
return( (x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2) );
}

void getplaquetteattribs( struct plaquette *Plaq, float u, float v, float w){
float del1x,del1y,del1z,del2x,del2y,del2z,dotprod,nlength,viewlength;
del1x=Plaq->vertex[1][0]-Plaq->vertex[0][0];
del1y=Plaq->vertex[1][1]-Plaq->vertex[0][1];
del1z=Plaq->vertex[1][2]-Plaq->vertex[0][2];
del2x=Plaq->vertex[3][0]-Plaq->vertex[0][0];
del2y=Plaq->vertex[3][1]-Plaq->vertex[0][1];
del2z=Plaq->vertex[3][2]-Plaq->vertex[0][2];
Plaq->norml[0]=del1y*del2z-del2y*del1z;
Plaq->norml[1]=del1z*del2x-del2z*del1x;
Plaq->norml[2]=del1x*del2y-del2x*del1y;
nlength=(Plaq->norml[0])*(Plaq->norml[0])+(Plaq->norml[1])*(Plaq->norml[1])+(Plaq->norml[2])*(Plaq->norml[2]);
nlength=sqrt(nlength);
viewlength=sqrt(u*u+v*v+w*w);
if(nlength!=0.0)
Plaq->face=((Plaq->norml[0])*u+(Plaq->norml[1])*v+(Plaq->norml[2])*w)/(nlength*viewlength);
else
Plaq->face=0.0;
Plaq->clr=(int)(65535.0*(1.0+tanh(2.0*Plaq->face))/2.0);
}


float function(float x, float y){
    /* this is the radiation function itself */
float f,factor,beta,KayEll;
 beta=0.6; /* speed of particle */
 factor=1.0-beta*cos(y); /* a universally used factor */

    /* linear acceleration */
 f=pow(sin(y),3.0)/pow(factor,5.0);

/* synchrotron radiation */
  f=(factor*factor-(1.0-beta*beta)*sin(y)*sin(y)*cos(x)*cos(x));
    f=f*f/pow(factor, 5.0); 

/* antennas */

/* center-fed */
  KayEll=5.0;
    f=(cos(KayEll*cos(y))-cos(KayEll))/sin(y); 
    f=f*f; 

/* full wave */
 /* f=sin(PI*cos(y))/sin(y);
    f=f*f; */
return(f);
}

 int mod(int A, int B){
     /* quick and dirty fix to impose P.B.C's */
  if(A>=B)
return(A-B);
  else
 return(A-B*(A/B));
  }