#include<stdio.h>
#include<string.h>
#include<math.h>
#include<xygraphm.h>

#define TITLE   ("Ray paths for the IGRF model")
//#define TITLE   ("Ray paths for the dipole field model")

#define MAX_RAY (20)
#define NLPMAX  (30000)
#define RAD     (M_PI/180)
#define RE      (6370.0)

int main(){

  int i,j,ipage;
  double tmp,tmp2;

  int nray;            // number of rays
  int istep[MAX_RAY];  // number of step of each ray

  FILE *fp;
  int str_num;
  char str[128];
  double thi[MAX_RAY][NLPMAX], rrr[MAX_RAY][NLPMAX], phi[MAX_RAY][NLPMAX];
  double del[MAX_RAY][NLPMAX], eps[MAX_RAY][NLPMAX], dns[MAX_RAY][NLPMAX];
  double ref[MAX_RAY][NLPMAX], dtm[MAX_RAY][NLPMAX];
  double ph0[MAX_RAY][NLPMAX], ph1[MAX_RAY][NLPMAX];
  
  double ex[181], ey[181];          // earth
  double qx[2], qy[2];              // equator
  double fx[4][181], fy[4][181];    // geomagnetic field lines
  double rx[NLPMAX], ry[NLPMAX];    // ray path
  double kk[NLPMAX];                // k-vector

  // XYGRAPH
  XYGRAPH xyg;
  double xf, yf, xal, yal;
  int kc,mode;
  double xmin, xmax, ymin, ymax, h, g, gt, gs;
  char tx[200], ty[200], ttl[200], stl[200];
  int itx, ity, ittl, istl;
  char comment[200];
  sprintf(comment,"%s",TITLE);


  // Data read
  if((fp=fopen("fort.21","r"))!=NULL){
    nray = 0;
    for(i=0; i<MAX_RAY; i++) istep[i]=0;
    while( fgets(str, sizeof(str), fp) != NULL ){   
      str_num = sscanf(str,"%lf %lf %lf %lf %lf %lf %lf %lf",
		       thi[nray]+istep[nray], rrr[nray]+istep[nray], phi[nray]+istep[nray], 
		       del[nray]+istep[nray], eps[nray]+istep[nray], dns[nray]+istep[nray], 
		       ref[nray]+istep[nray], dtm[nray]+istep[nray]);
      if (str_num!=8) {
	if (istep[nray]<2) istep[nray]=0; else nray++;
      }
      istep[nray]++;
    }
    fclose(fp);
  } else {
    fprintf(stderr,"fort.21 cannot be opened\n");
    return -1;
  }
  for(i=0; i<MAX_RAY; i++) istep[i]--;



  for(ipage=0; ipage<nray; ipage++){

    xyg.baseh->MSymbl(5.0, 24.5, 0.9, comment, 0.0, strlen(comment));

    // Axis, Label, Title 
    xf   =  3.0;
    yf   =  3.0;
    xal  = 20.0;
    yal  = 20.0;
    kc   =  5;
    xmin =  0.0;
    xmax =  4.0;
    ymin = -2.0;
    ymax =  2.0;
    h    =  0.5;
    g    =  0.5;
    sprintf(tx,"[Re]");
    sprintf(ty,"[Re]");

    xyg.Colpen(1);
    xyg.RnwPen(1);
    xyg.XYAx(xf, yf, xal, yal, kc, xmin, xmax, ymin, ymax, h,
	     tx, strlen(tx), ty, strlen(ty), g);

    // Earth 
    for(i=0; i<181; i++){
      tmp = (double)(i-90)*RAD;
      ex[i]=cos(tmp);
      ey[i]=sin(tmp);
    }
    xyg.XYLine(ex,ey,181);

    // Latitude
    mode = 0;
    xyg.baseh->DPlot(qx[0],qy[0],mode);
    for(j=-80; j<=80; j+=10){
      tmp = (double)(j)*RAD;
      qx[0] = cos(tmp);
      qy[0] = sin(tmp);
      qx[1] = 6.0*cos(tmp);
      qy[1] = 6.0*sin(tmp);
      xyg.XYLine(qx,qy,2);
    }
      
    // Geomagnetic field lines for L=5,4,3,2
    for(i=0; i<4; i++){
      for(j=0; j<=180; j++){
	tmp = cos((double)(j-90)*RAD);
	if((double)(i+2)*tmp*tmp>=1.0){
	  tmp2 = cos((double)(j-90)*RAD);
	  fx[i][j]=(i+2)*tmp2*tmp2*cos((double)(j-90)*RAD);
	  fy[i][j]=(i+2)*tmp2*tmp2*sin((double)(j-90)*RAD);
	}else{
	  fx[i][j]=-1.0;
	  fx[i][j]=-1.0;
	}
      }
    }
    xyg.XYMiss((int)(-1.0));
    for(i=0; i<4; i++) xyg.XYLine(fx[i],fy[i],181);
    mode = 1;
    xyg.baseh->DPlot(fx[0][0],fy[0][0],mode);

    // Ray path & k-vector
    for(i=0; i<istep[ipage]; i++){
      rrr[ipage][i+1] = rrr[ipage][i+1] + RE;  // Altitude -> Radial distance
      rx[i] = rrr[ipage][i+1]/RE*cos(thi[ipage][i+1]*RAD);
      ry[i] = rrr[ipage][i+1]/RE*sin(thi[ipage][i+1]*RAD);
      kk[i] = thi[ipage][i+1]-del[ipage][i+1];
    }
    xyg.Colpen(6);
    xyg.RnwPen(6);
    xyg.XYLine(rx, ry, istep[ipage], 15, 0.1, 1);
    
#define DECIM  (20)
    xyg.Colpen(2);
    xyg.RnwPen(2);
    for(i=0; i<(int)(istep[ipage]/DECIM); i++){
      xyg.ksl->XYArow(rx[i*DECIM], ry[i*DECIM], 0.5, kk[i*DECIM]);
    }

    // Electron density
    xf   = 26.0;
    yf   = 14.5;
    xal  = 10.0;
    yal  =  8.5;
    kc   = 2005;
    xmin =  0.0;
    xmax =  0.1;
    ymin =1.0e2;
    ymax =1.0e6;
    h    =  0.4;
    g    =  0.4;
    sprintf(tx,"Time [sec]");
    sprintf(ty,"Electron density [/cm|3]");

    xyg.Colpen(1);
    xyg.RnwPen(1);
    xyg.XYAx(xf, yf, xal, yal, kc, xmin, xmax, ymin, ymax, h,
	     tx, strlen(tx), ty, strlen(ty), g);
    for(i=0; i<istep[ipage]; i++){
      dns[ipage][i] = dns[ipage][i]/8.98;
      dns[ipage][i] = dns[ipage][i]*dns[ipage][i];
    }
    xyg.XYLine(dtm[ipage], dns[ipage], istep[ipage]);

    // Refractive index
    xf    = 26.0;
    yf    =  3.0;
    xal   = 10.0;
    yal   =  8.5;
    kc    = 2005;
    xmin  = 134.5;
    xmax  = 135.5;
    ymin  = 1.0e0;
    ymax  = 1.0e4;
    h     = 0.4;
    g     = 0.4;
    sprintf(tx,"Time [sec]");
    sprintf(ty,"Refractive index");
    xyg.Colpen(1);
    xyg.RnwPen(1);
    xyg.XYAx(xf, yf, xal, yal, kc, xmin, xmax, ymin, ymax, h,
	     tx, strlen(tx), ty, strlen(ty), g);
    xyg.XYLine(dtm[ipage], ref[ipage], istep[ipage]);

    xyg.NewPage();
  }

  return 0;
}