#include </usr/local/pgplot/cpgplot.h>
/* change from </usr/local/pgplot/cpgplot.h> to the heada file in your system*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifndef EXIT_FAILURE
#define EXIT_FAILURE 1
#endif
#ifndef EXIT_SUCCESS
#define EXIT_SUCCESS 0
#endif
void plot1();
main()
 {
  char y;
  if(cpgbeg(0, "/ps", 1, 10) != 1)
  exit(EXIT_FAILURE);
  cpgask(1);
  plot1();
  cpgend();
  printf("PRINT ok?(y/n)\n");
  scanf("%s",&y);
   if(y=='y'){
  system("lp pgplot.ps");}
 /* system("lp pgplot.ps");*/
 printf("END\n");
  return EXIT_SUCCESS;
}
void plot1()
{
  int  REV,Y,M,D,a=0,j,k=0,i,y,m,d=1,t=0,l=0,kk=0,kkk=0,ihh,imm,mm;
  int p,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15,s16,
      s17,s18,s20,s21,s22,dem,t4[5500],n4[5500],L,Te,Ne,TE[13][48],NE[13][48];
  long s19;
  char fname1[50],fname2[50],fname3[50],al,buf[250],MO[7],FL[5];
  double dx[13][48],dy[13][48],dxn[13][48],dyn[13][48],DX,DY,MLAT,T,LT,DXn,DYn,
         t0[5500],n0[5500],t1[5500],t2[5500],t3[5500],n1[5500],n2[5500],n3[5500],
         t5[5500],t6[5500],t7[5500],t8[5500];
  float xp[2],yp[2],xr[5500], yr[5500],yrr[5500];
  float xmin,xmax,ymin,ymax,xd,x1,y1,llx[2],lly[2];
  char  string[10],string1[10],string2[10],string3[10],string4[10],
	string5[10],string6[10];
   FILE *fa,*fb,*fc,*fd;
   
   while(l<=12){while(t<=47){TE[l][t]=0;dx[l][t]=0.0;dy[l][t]=0.0;NE[l][t]=0;dxn[l][t]=0.0;dyn[l][t]=0.0;t++;}t=0;l++;}
  t=0;l=0;
 
  printf("input Data file month. ex 8102\n");
  scanf("%s",&FL);
  mm=atoi(FL);
  if((mm<=8112 && mm>=8111 ) || mm==8201){strcpy(MO,"1112-1");/*MO='1112-1';*/}
  if((mm<=8104 && mm>=8102) || (mm<=8204 && mm>=8202) ){strcpy(MO,"234");}
  if((mm<=8107 && mm>=8105) || (mm<=8207 && mm>=8205) ){strcpy(MO,"567");}
  if(mm<=8110 && mm>=8108 ){strcpy(MO,"8910");}
 
 printf("Which Average file? High solar activity(F10.7>175)=1 or Low solar activity(F10.7<200)=2.\n");
  scanf("%d",&p);
   if(p==2){sprintf(fname1,"%s%s","./AVE/",MO);}
   if(p==1){ sprintf(fname1,"%s%s%s","./AVE/",MO,"h");}
   if((fb=fopen(fname1,"r"))==NULL){printf("TE ave AKIMASEN\n");exit(1);}
   printf("Te ave ok\n");
    while (fscanf(fb,"%d %lf %d %lf %lf",&L,&T,&Te,&DX,&DY)!=EOF){
    while(l<=12){
	    while(t<=47){
		if(L==-30+5*l && T==0.5*t){
			TE[l][t]=Te;
			dx[l][t]=DX;
			dy[l][t]=DY;
			}
			t++;
			}
		t=0;
		l++;
	}l=0;	}
	t=0;
	l=0;
  fclose(fb);
  if(p==2){sprintf(fname2,"%s%s","./AVE/n",MO);}
  if(p==1){sprintf(fname2,"%s%s%s","./AVE/n",MO,"h");}
   if((fd=fopen(fname2,"r"))==NULL){printf("NE ave AKIMASEN\n");exit(1);}
  printf("Ne ave ok\n");
    while (fscanf(fd,"%d %lf %d %lf %lf",&L,&T,&Ne,&DXn,&DYn)!=EOF){
    while(l<=12){
	    while(t<=47){
		if(L==-30+5*l && T==0.5*t){
			NE[l][t]=Ne;
			dxn[l][t]=DXn;
			dyn[l][t]=DYn;}
			t++;
			}
		t=0;
		l++;
	}l=0;	}
	t=0;
	l=0;
  fclose(fd);
 /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
   /*change the directry */
   sprintf(fname3,"%s%s","/Data/Data/Hinotori/hinotori.",FL);

 /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
   if ((fa=fopen(fname3,"r"))==NULL){printf("AKIMASEN\n");exit(1);} printf(" ok\n");
   printf("Orbit number?\n");
   scanf("%d",&REV);
 
 while(a<5499){t0[a]=t1[a]=t2[a]=t3[a]=n0[a]=n1[a]=n2[a]=n3[a]=0.0;t5[a]=t7[a]=t8[a]=0;a++;}
 
 /*Read Data and make averages Te*/
 while((fscanf(fa,"%d %d %d %d %d %d %d %d %d %d
		   %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d",
  &s1,&s2,&s3,&s4,&s5,&s6,&s7,&s8,&s9,&s10,&s11,&s12,&s13,&s14,&s16,&s15,&s17,
  &s18,&s19,&s20,&s21,&s22,&dem,&dem,&dem,&dem,&dem,&dem,&dem,&dem))!=EOF){
 if(s1==REV && s15<3250 && s15>=-3250)
 {if(k==0){ Y=s2; M=s3; D=s4;}
 if(s21!=0){ t=0;l=0;
  while(l<=12){
     while(t<=47){
	if((s15<-3000+500*l+250 && s15>=-3000+500*l) && (s11*60+s12<30*t+15 && s11*60+s12>=30*t )){
	t0[k]=s5+s6/60.0;
	t1[k]=s15/100.0;
	t2[k]=(s11*60+s12)/60.0;
	t3[k]=TE[l][t]+dx[l][t]*(s15/100.0-(-30.0+5.0*l))+dy[l][t]*((s11*60.0+s12)-30.0*t);
	t4[k]=s21;t5[k]=s14/100.0;t7[k]=s16/100.0;t8[k]=s17/100.0;k++;}
				
	if((s15<-3000+500*l+250 && s15>=-3000+500*l) &&(s11*60+s12<30*t && s11*60+s12>=30*t-15) ){/*printf("%d\n",k);*/
	t0[k]=s5+s6/60.0;
	t1[k]=s15/100.0;
	t2[k]=(s11*60+s12)/60.0;
	t3[k]=TE[l][t]+dx[l][t]*(s15/100.0-(-30.0+5.0*l))-dy[l][t-1]*(30.0*t-(s11*60.0+s12));
	t4[k]=s21;t5[k]=s14/100.0;t7[k]=s16/100.0;t8[k]=s17/100.0;k++;}
				 
	if((s15<-3000+500*l && s15>=-3000+500*l-250) &&(s11*60+s12<30*t+15 && s11*60+s12>=30*t) ){/*printf("%d\n",k);*/
	t0[k]=s5+s6/60.0;
	t1[k]=s15/100.0;
	t2[k]=(s11*60+s12)/60.0;
	t3[k]=TE[l][t]-dx[l-1][t]*((-30.0+5.0*l)-s15/100.0)+dy[l][t]*((s11*60.0+s12)-30.0*t);
	t4[k]=s21;t5[k]=s14/100.0;t7[k]=s16/100.0;t8[k]=s17/100.0;k++;}
				 
	if((s15<-3000+500*l && s15>=-3000+500*l-250) &&(s11*60+s12<30*t && s11*60+s12>=30*t-15 )){/*printf("%d\n",k);*/
	t0[k]=s5+s6/60.0;
	t1[k]=s15/100.0;
	t2[k]=(s11*60+s12)/60.0;
	t3[k]=TE[l][t]-dx[l-1][t]*((-30.0+5.0*l)-s15/100.0)-dy[l][t-1]*(30.0*t-(s11*60.0+s12));
	t4[k]=s21;t5[k]=s14/100.0;t7[k]=s16/100.0;t8[k]=s17/100.0;k++;}			 
	t++;}
	t=0;
    l++;
    } l=0;

}
  if(s19!=0){t=0;l=0;
  while(l<=12){
	while(t<=47){
	if((s15<-3000+500*l+250 && s15>=-3000+500*l) && (s11*60+s12<30*t+15 && s11*60+s12>=30*t )){
	n0[kk]=s5+s6/60.0;
	n1[kk]=s15/100.0;
	n2[kk]=(s11*60+s12)/60.0;
	n3[kk]=NE[l][t]+dxn[l][t]*(s15/100.0-(-30.0+5.0*l))+dyn[l][t]*((s11*60.0+s12)-30.0*t);n4[kk]=s19;kk++;}
				
	if((s15<-3000+500*l+250 && s15>=-3000+500*l) &&(s11*60+s12<30*t && s11*60+s12>=30*t-15) ){
	n0[kk]=s5+s6/60.0;
	n1[kk]=s15/100.0;
	n2[kk]=(s11*60+s12)/60.0;
	n3[kk]=NE[l][t]+dxn[l][t]*(s15/100.0-(-30.0+5.0*l))-dyn[l][t-1]*(30.0*t-(s11*60.0+s12));n4[kk]=s19;kk++;}
				 
	if((s15<-3000+500*l && s15>=-3000+500*l-250) &&(s11*60+s12<30*t+15 && s11*60+s12>=30*t) ){
	n0[kk]=s5+s6/60.0;
	n1[kk]=s15/100.0;
	n2[kk]=(s11*60+s12)/60.0;
	n3[kk]=NE[l][t]-dxn[l-1][t]*((-30.0+5.0*l)-s15/100.0)+dyn[l][t]*((s11*60.0+s12)-30.0*t);n4[kk]=s19;kk++;}
				 
	if((s15<-3000+500*l && s15>=-3000+500*l-250) &&(s11*60+s12<30*t && s11*60+s12>=30*t-15 )){
	n0[kk]=s5+s6/60.0;
	n1[kk]=s15/100.0;
	n2[kk]=(s11*60+s12)/60.0;
	n3[kk]=NE[l][t]-dxn[l-1][t]*((-30.0+5.0*l)-s15/100.0)-dyn[l][t-1]*(30.0*t-(s11*60.0+s12));n4[kk]=s19;kk++;}			 
	t++;}
	t=0;
    l++;
    }
    l=0; 
    }

 } l=0;t=0;
 if(s1>REV){break;} } 
  printf("%d %f\n",k,t0[k]);
  fclose(fa);
 /* set */
 xmin = 0.0;
 xmax = 24.0;
  cpgsch(8.0);
 
/* 1 */
  ymin = 4.0;
  ymax = 7.0;
  cpgenv(xmin,xmax,ymin,ymax, 0, -2);
  cpgsvp(0.2,1.0,-2.0,0.4);
  if(p==2){ sprintf(string4,"%d %d %d Orbit Number %d",Y,M,D,REV);}
  if(p==1){ sprintf(string4,"%d %d %d Orbit Number %d high S.A",Y,M,D,REV);}
  cpglab(" ", "(Ne)",string4 );
  for(i=0; i<kk; i++){
    yrr[i]= log10(n3[i]);
    yr[i]= log10(n4[i]);
    xr[i] = n2[i];
    
  }
  xp[0]=0.0; xp[1]=24.0;
  yp[0]=0.0; yp[1]=0.0;
  cpgpt(kk, xr, yr,5);
  cpgpt(kk, xr, yrr,1);
  cpgline(2,xp,yp);
  cpgbox("BCTS",0.0,0,"BCTSNV",0.0,0);
/* 2 */
  ymin = 800.0;
  ymax = 4200.0;
  cpgenv(xmin,xmax,ymin,ymax, 0, -2);
  cpgsvp(0.2,1.0,-3.7,-1.3);
  cpglab(" ", "(Te)",""  );
  for(i=0; i<k; i++) {
   yr[i] =t4[i];
   yrr[i]= t3[i];
   xr[i] = t2[i];
  }
  cpgpt(k, xr,yr,5);
  cpgpt(k, xr, yrr,1);
  cpgbox("BCTS",0.0,0,"BCTSNV",0.0,0);
  llx[0]=0.0;llx[1]=24.0;lly[0]=lly[1]=0.0; 
/* 3 */
 ymin = -27.5;
  ymax = 27.5;
  cpgenv(xmin,xmax,ymin,ymax, 0, -2);
  cpgsvp(0.2,1.0,-4.5,-3.0);
  cpglab(" ", "MLAT", " ");
   for(i=0; i<k; i++) {
     yr[i] = t1[i];
    xr[i] = t2[i]; 
  }
  cpgpt(k, xr,yr,20);
  cpgline(2,llx,lly);
  cpgbox("BCTS",0.0,0,"BCTSNV",0.0,0);

/* parameters */
  ymin = 0.0;
  ymax = 10.0;
  cpgswin(xmin,xmax,ymin,ymax);
  cpgsvp(0.2,1.0,-7.0,-5.0);
  x1 = xmin - (xmax-xmin)/11.0;
  y1 = 10.0;
  cpgptxt(x1,y1,0.0,0.0,"UT");
  cpgptxt(x1,y1-1.5,0.0,0.0,"LT");
  cpgptxt(x1,y1-3.0,0.0,0.0,"LAT");
  cpgptxt(x1,y1-4.5,0.0,0.0,"MLAT");
  cpgptxt(x1,y1-6.0,0.0,0.0,"LONG");
  cpgptxt(x1,y1-7.5,0.0,0.0,"ALT");
 for(xd=xmin;xd<=xmax;xd=xd+2.0){
	 ihh = (int)(xd);
	 imm = (int)( (xd-ihh)*60.0 );
	 x1 = (float)ihh + (float)imm/60.0 - (xmax-xmin)/100.0;
         sprintf(string,"%2d:%2d\0",ihh,imm);
         cpgptxt(x1,y1-1.5,0.0,0.0,string);
         
 for(i=0; i<k; i++) {
 if(xd==(t2[i]) ){
 		 sprintf(string2,"%2d:%d\0",(int)(t0[i]),(int)((t0[i]-(int)(t0[i]))*60.0));
 		sprintf(string3,"%2.2lf\0",t7[i]);
 		sprintf(string1,"%2.2lf\0",t1[i]);
   		sprintf(string5,"%3.2lf\0",t5[i]);
		sprintf(string6,"%3.2lf\0",t8[i]);
                cpgptxt(x1,y1,0.0,0.0,string2);
 		cpgptxt(x1,y1-3.0,0.0,0.0,string3);
 		cpgptxt(x1,y1-4.5,0.0,0.0,string1);
		cpgptxt(x1,y1-6.0,0.0,0.0,string5);
		cpgptxt(x1,y1-7.5,0.0,0.0,string6);
 		break;}

 
 if(xd!=(t2[i]) && (xd<=t2[i]+5.0/60.0 && xd>=t2[i]-5.0/60.0)){
 		sprintf(string2,"%2d:%d\0",(int)(t0[i]),(int)((t0[i]-(int)(t0[i]))*60.0));
		sprintf(string3,"%2.2lf\0",t7[i]);
 		sprintf(string1,"%2.2lf\0",t1[i]);
		sprintf(string5,"%3.2lf\0",t5[i]);
		sprintf(string6,"%3.2lf\0",t8[i]);
               	cpgptxt(x1,y1,0.0,0.0,string2);
 		cpgptxt(x1,y1-3,0.0,0.0,string3);
 		cpgptxt(x1,y1-4.5,0.0,0.0,string1);
		cpgptxt(x1,y1-6.0,0.0,0.0,string5);
		cpgptxt(x1,y1-7.5,0.0,0.0,string6);
		break;}
if(i==k-1) {	cpgptxt(x1,y1,0.0,0.0,"***");
 		cpgptxt(x1,y1-3,0.0,0.0,"***");		
 		cpgptxt(x1,y1-4.5,0.0,0.0,"***");
 		cpgptxt(x1,y1-6.0,0.0,0.0,"***");
 		cpgptxt(x1,y1-7.5,0.0,0.0,"***");}
		}
  }
  return;
}