/* { make ; pag; ash /users/pbl/BRCO04[100:400] +y 10000 -ly } */
/* This program draws a spectrum, possibly with error bars
   possible parameters:
   ash sp[1:1024] {-y 0} {+y 1000} {-b} {-v} {-c}
       sp[min:max] specifies spectrum name and range (default is stdin)
       -y %E specifies the lower bound
       +y %E specifies the upper bound
       -lx logarithmic x-scale
       -ly logarithmic y-scale
       -d display only (don't draw axis)
       -pag clear screen before drawing
       -use use old parameters
       -cur prints cursor position after drawing
       -b draws error bars \
       -v draws vectors     - default is point only
       -c draws circles    /
   The boundaries of the real data area of the screen are saved
   to /tmp/ashbdry as ASCII integers

                                          Programmer: RAKO
*/

#include <stdio.h>
#ifdef AMIGA
#include <gfxamiga.h>
#include <workbench/workbench.h>
#include <workbench/startup.h>
#include <workbench/icon.h>
extern struct WBStartup *WBenchMsg;
extern struct IconBase *IconBase;
#else
#include <gfx.h>
#endif

int   dot_flg=TRUE,
      bar_flg=FALSE,
      vec_flg=FALSE,
      point_style=0,
      point_size=0,
      auto_flg=TRUE,
      logx=FALSE,
      logy=FALSE,
      vecpos=0;
float *spc,*err,*tim;

help()
{
printf("Draw a spectrum with axis\n");
printf(" possible parameters:\n");
printf(" ash sp[1:1024] {-y 0} {+y 1000} {-b} {-v} {-c}\n");
printf("  sp[min:max] specifies spectrum name and range (default is stdin)\n");
printf("  -y n.m    specifies the lower bound\n");
printf("  +y n.m    specifies the upper bound\n");
printf("  -lx       logarithmic x-scale\n");
printf("  -ly       logarithmic y-scale\n");
printf("  -d        display only (don't draw axis)\n");
printf("  -pag      clear screen before drawing\n");
printf("  -use      use old parameters\n");
printf("  -cur      prints cursor position after drawing\n");
printf("  -b        draws error bars \\\n");
printf("  -v        draws vectors     - default is point only\n");
printf("  -c        draws circles    /\n");
printf("  -cir n    draws circles with radius n.m\n");
printf("  -crf n    draws filled circles with radius n\n");
printf("  -tri n    draws triangles with baseline n (negative n gives\n");
printf("                                            inverted triangles)\n");
printf("  -trf n    draws filled triangles with baseline n\n");
printf("  -qua n    draws quadrates with baseline n\n");
printf("  -qaf n    draws filled quadrates with baseline n\n");
printf("  -mb n     defines bottom margin\n");
printf("  -mt n     defines top margin\n");
printf("  -mr n     defines right margin\n");
printf("  -ml n     defines left margin\n");
}

minimax(spc,xmax,min,max)
float spc[],*min,*max;
int   xmax;
{  /* findig maximum and minimum of spectrum */
int n;
   *max = -1.0E10; *min = 1.0E10;
   for(n=0;n<xmax;n++) {
      if(spc[n] > *max) *max=spc[n];
      if(spc[n] < *min) *min=spc[n];
   }
   if(*max == *min) {
      *min = 0.0;
      *max = 10.0;
   }
   if(*max == -1.0E10) *max = 10.0;
   if(*min == 1.0E10) *min = 0.0;
}

disp(y,e,t)       /* plot a point in the window */
float y,e,t;
{
int   n,m,xx,yy,ee;
float fx1,fx2,fy1,fy2,fr,phi,finc;

   xx=t; yy=y; ee=e;
   if(dot_flg) {
      posita(xx,yy); vectoa(xx,yy);
   }
   if(bar_flg) {
      posita(xx,yy-ee); vectoa(xx,yy+ee);
      posita(xx,yy);
   }
   if(vec_flg) {
      if(vecpos==0) {
         vecpos=1;
         posita(xx,yy);
      }
      vectoa(xx,yy);
   }
   switch(point_style) {
   case 0:
      return(0);
      break;
   case 1:     /* circle */
      fr=point_size;
      xx = t ; yy = y+fr; posita(xx,yy);
      finc=1.57079/fr;
      phi=finc;
      fx1=0; fy1=fr;
      for(n=0; n < (4*point_size); n++) {
         fx2=fr*sin(phi); fy2=fr*cos(phi);
         xx = t + fx2; yy = y + fy2; vectoa(xx,yy);
         phi=phi+finc;
      }
      xx=t; yy=y; posita(xx,yy);
      break;
   case 2:     /* filled circle */
      fr=point_size;
      phi=0.0;
      finc=1.57079/fr;
      fx1=0; fy1=fr;
      for(n=0;n <= (point_size + 1);n++) {
         fx2=fr*sin(phi); fy2=fr*cos(phi);
         xx=(t+fx1); yy=(y+fy1); posita(xx,yy);
         xx=(t-fx1); vectoa(xx,yy);
         yy=(y-fy1); posita(xx,yy);
         xx=(t+fx1); vectoa(xx,yy);
         fx1=fx2; fy1=fy2;
         phi=phi+finc;
      }
      xx=t; yy=y; posita(xx,yy);
      break;
   case 3:     /* triangle */
      xx=t+point_size; yy=y-point_size; posita(xx,yy);
      xx=t-point_size; vectoa(xx,yy);
      yy=y+point_size; xx=t; vectoa(xx,yy);
      xx=t+point_size; yy=y-point_size; vectoa(xx,yy);
      xx=t; yy=y; posita(xx,yy);
      break;
   case 4:     /* filled triangle */
      xx=t-point_size;
      for(n=0;n<(2*point_size);n++) {
         yy=y+point_size; posita((int)t,yy);
         yy=y-point_size; vectoa(xx,yy);
         xx=xx+1;
      }
      xx=t; yy=y; posita(xx,yy);
      break;
   case 5:     /* quadrate */
      xx=t+point_size; yy=y-point_size; posita(xx,yy);
      xx=t-point_size; vectoa(xx,yy);
      yy=y+point_size; vectoa(xx,yy);
      xx=t+point_size; vectoa(xx,yy);
      yy=y-point_size; vectoa(xx,yy);
      xx=t; yy=y; posita(xx,yy);
      break;
   case 6:     /* filled quadrate */
      yy=y-point_size;
      for(n=0;n<(2*point_size);n++) {
         xx=t+point_size; posita(xx,yy);
         xx=t-point_size; vectoa(xx,yy);
         yy=yy+1;
      }
      xx=t; yy=y; posita(xx,yy);
      break;
   }
}

float slog(x)
float x;
{
float erg,y;
erg = 0.0;
if( x > 0.0) erg = log(x);  /* / 2.3026; */
return(erg);
}

float pow10intlog10(x) /* same as n=pow(10,floor(log10( (double)d + 1.0E-10))) */
float x;
{
int ix,n;
float xx;

   xx=1.0;
   if(x > 1.0) {
      while(xx < x) xx = xx * 10.0;
      xx = xx / 10.0;
   }
   if(x < 1.0) {
      while(xx > x) xx = xx / 10.0;
   }
   return(xx);
}


main(argc,argv)
int argc;
char *argv[];
{
int   n,m,i,j,xmax;
float x,y,
      ymax = 1.0,
      ymin = 0.0,
      yfak = 1.0,
      tmin = 0.0,
      tmax = _MAXSPCLEN,
      tfak = 1.0,
      xoffset = 0,
      yoffset = 0;
int   leftmargin = 100,
      bottommargin = 30,
      rightmargin = _TEKXMAX-30,
      topmargin = _TEKYMAX-30;
char  z[80],comment[80];
FILE  *fp;
#ifdef AMIGA
struct DiskObject *Lock;
struct WBArg *Arg;
#endif

   if(checkopt(argc,argv,"-lx",z)) logx=TRUE;
   if(checkopt(argc,argv,"-ly",z)) logy=TRUE;

   if(checkopt(argc,argv,"-y",z)) {
      auto_flg=FALSE;
      ymin=atosf(z);
   }
   if(checkopt(argc,argv,"+y",z)) {
      auto_flg=FALSE;
      ymax=atosf(z);
   }
   if(checkopt(argc,argv,"-b",z)) {
      dot_flg=FALSE;
      bar_flg=TRUE;
   }
   if(checkopt(argc,argv,"-v",z)) {
      dot_flg=FALSE;
      vec_flg=TRUE;
   }
   if(checkopt(argc,argv,"-c",z)) {
      dot_flg=FALSE;
      point_style=1;
      point_size=4;
   }
   if(checkopt(argc,argv,"-cir",z)) {
      dot_flg=FALSE;
      point_style=1;
      point_size=atoi(z);
   }
   if(checkopt(argc,argv,"-crf",z)) {
      dot_flg=FALSE;
      point_style=2;
      point_size=atoi(z);
   }
   if(checkopt(argc,argv,"-tri",z)) {
      dot_flg=FALSE;
      point_style=3;
      point_size=atoi(z);
   }
   if(checkopt(argc,argv,"-trf",z)) {
      dot_flg=FALSE;
      point_style=4;
      point_size=atoi(z);
   }
   if(checkopt(argc,argv,"-qua",z)) {
      dot_flg=FALSE;
      point_style=5;
      point_size=atoi(z);
   }
   if(checkopt(argc,argv,"-qaf",z)) {
      dot_flg=FALSE;
      point_style=6;
      point_size=atoi(z);
   }
   if(checkopt(argc,argv,"-mt",z)) topmargin = atoi(z);
   if(checkopt(argc,argv,"-mb",z)) bottommargin = atoi(z);
   if(checkopt(argc,argv,"-ml",z)) leftmargin = atoi(z);
   if(checkopt(argc,argv,"-mr",z)) rightmargin = atoi(z);
   
   tekopen();

   if(checkopt(argc,argv,"-pag",z)) dperas();

   spc= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   err= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   tim= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   if(tim==NULL) {
      printf("sorry, not enough memory\n");
      exit(-1);
   }
   if(argc > 1) strcpy(z,argv[1]);

#ifdef AMIGA
   if(argc == 0) {
      IconBase = (struct IconBase *) OpenLibrary("icon.library", 0L);
      Arg = WBenchMsg->sm_ArgList; Arg++;
      CurrentDir(Arg->wa_Lock);
      strcpy(z,Arg->wa_Name);
      if(strlen(z) < 2) exit(0);
      dperas();
   }
#endif

   xmax=readspec(z,spc,err,tim,comment);
   if(auto_flg) {       /* findig maximum and minimum of spectrum */
      minimax(spc,xmax,&ymin,&ymax);
      ymax = 1.1 * ymax;
      if(ymin < 0.0)  {
         ymin = 1.1 * ymin;
      } else {
         ymin = 0.9 * ymin;
      }
   }
   minimax(tim,xmax,&tmin,&tmax);


   if(logy) {
      for(n=0;n<xmax;n++) {
         spc[n] = slog(spc[n]);
         err[n] = slog(err[n]);
      }
   ymin = slog(ymin);
   ymax = slog(ymax);
   }

   if(logx) {
      for(n=0;n<xmax;n++) {
         tim[n] = slog(tim[n]);
      }
   tmin = slog(tmin);
   tmax = slog(tmax);
   }



   if(checkopt(argc,argv,"-use",z)) {
printf("ymax=%E\n",ymax);
      fp=fopen(ASHBDRY,"r");
      fscanf(fp,"%d\n%d\n%d\n",&leftmargin,&bottommargin,&rightmargin);
      fscanf(fp,"%d\n%d\n%E\n",&topmargin,&_win_flg,&tmin);
      fscanf(fp,"%E\n%E\n%E\n",&tmax,&ymin,&ymax);
printf("ymax=%E\n",ymax);
      fclose(fp);
   }

   yfak=(topmargin-bottommargin)/(ymax-ymin);
   tfak=(rightmargin-leftmargin)/(tmax-tmin);
   xoffset= (- tmin) *tfak + leftmargin;
   yoffset= (- ymin) *yfak + bottommargin;

   if(!checkopt(argc,argv,"-d",z)) {
      axis(leftmargin,bottommargin,rightmargin,topmargin,tmin,tmax,ymin,ymax);
   }


   for(n=0;n<xmax;n++) {
      disp(spc[n] * yfak + yoffset , err[n] * yfak , tim[n] * tfak + xoffset);
   }

   if(!checkopt(argc,argv,"-d",z)) {
      posita(leftmargin+20,topmargin+8);
      gfxtext(comment,0.0);
      gfxflush();
   }

   unlink(ASHBDRY);
   /* print out boundaries of usable data area */
   fp=fopen(ASHBDRY,"w");
   fprintf(fp,"%d\n%d\n%d\n",leftmargin,bottommargin,rightmargin);
   fprintf(fp,"%d\n%d\n%E\n",topmargin,_win_flg,tmin);
   fprintf(fp,"%E\n%E\n%E\n",tmax,ymin,ymax);
   fprintf(fp,"%E\n",_tica);
   fclose(fp);
   close(_tek4014);
   free(spc); free(err); free(tim);
   exit(0);
}

/* ----------------------------------------------------------------------
   Draw axis and numbers
   ---------------------------------------------------------------------- */

prettystr(s,x) /* generate string from fp number */
char s[];
float x;
{
int ix,fx;

   sprintf(s,"%-3.3f",x);
   if(fabs(x) < 0.001) sprintf(s,"%-.3E",x);
   if(fabs(x) > 9999) sprintf(s,"%-1.1E",x);
   ix=x; fx=ix; if((x - fx) == 0.0) sprintf(s,"%-d",ix);
   if(fabs(x)<1.0E-5) strcpy(s,"0");

}

axis(sx,sy,ex,ey,xmin,xmax,ymin,ymax) /* draw frame, ticks and numbers */
int sx,sy,ex,ey;
float xmin,xmax,ymin,ymax;
{
int n,i,x,y;
float nx,ny,ninc;
char s[80];

   /* draw frame */
   posita(sx,sy);
   vectoa(ex,sy);
   vectoa(ex,ey);
   vectoa(sx,ey);
   vectoa(sx,sy);


   x = sx; y = sy;
      
   calcnums(xmin,xmax,&nx,&ninc); x = -1000;
   if(logx) nx = pow10intlog10(xmin);
   while(x < ex) {
      x = sx + ((ex-sx) * ((nx-xmin)/(xmax-xmin)));
      if(logx) x = sx + ((ex-sx) * ((slog(nx)-xmin)/(xmax-xmin)));
      if(x>ex) break;
      posita(x,sy); vectoa(x,sy+5);
      posita(x,ey); vectoa(x,ey-5);
      prettystr(s,nx);
      posita(x-(textlen(s)/2),0);
      gfxtext(s,0.0);
      if(logx) {
         if(nx == 0) nx = 0.1;
         nx = nx * 10;
      } else {
         nx = nx + ninc;
      }
    }

   calcnums(ymin,ymax,&ny,&ninc); y = -1000;
   if(logy) ny = pow10intlog10(ymin);
   while(y < ey) {
      y = sy + ((ey-sy) * ((ny-ymin)/(ymax-ymin)));
      if(logy) y = sy + ((ey-sy) * ((slog(ny)-ymin)/(ymax-ymin)));
      if(y>ey) break;
      posita(sx,y); vectoa(sx+5,y);
      posita(ex,y); vectoa(ex-5,y);
      prettystr(s,ny);
      posita(sx-(2*textlen(s)),y-(4 * _gfxfontsz));
      gfxtext(s,0.0);
      if(logy) {
         if(ny == 0) ny = 0.1;
         ny = ny * 10;
      } else {
         ny = ny + ninc;
      }
   }
}


calcnums(min,max,ny,ninc) /* calculate increment for axis numbers */
float min,max,*ny,*ninc;
{
float d,n;
int m,y;

   d=fabs(max-min);
   n=pow10intlog10(d);
   m=d/n;
   switch(m) {
   case 1:
      *ninc=n/10;
      break;
   case 2:
      *ninc=n/4;
      break;
   default:
      *ninc=n;
      break;
   }
   *ny=n * floor(min/n);
   while(*ny < min) *ny = *ny + *ninc;
   if((d / *ninc) > 10.0) *ninc = *ninc * 2;
}