/*
   Subroutine: abc_gen.c

      This subroutine generates particle trajectories for the pre-defined
   generator (ig), and can be executed from the command-line or display
   window.  The routine stores particle track points in the trk*** arrays,
   and if a detector is entered, outputs hit data to standard output and fills
   the hit arrays hit*** (for display).
*/

#include "abc.h"

void abc_gen(int ig)
{
   int i,idt;
   int id,it,ip,is,io,i1,i2,i3,ir,ih;
   double mp,rm,rb,rl,pm,dp,dt,dx,dy,dz,x,y,z;
   double r1[4],r2[4],p1[4],p2[4];
   double r1d[4],r2d[4],p1d[4],p2d[4];

/* Check generator and track numbers */
   if(ig <= 0 || ig >= maxgen)
   {
      fprintf(stderr,"  Error: generator range (1,%d) \n",maxgen-1); return;
   }
   ip = genprt[ig]; if(ip <= 0)
   {
      fprintf(stderr,"  Error: undefined generator \n"); return;
   }
   if(trknum == maxtrk)
   {
      fprintf(stderr,"  Error: maximum track number (%d) \n",maxtrk); return;
   }

/* Initialize first step */
   r1[1] = genxpo[ig]; p1[1] = genxmo[ig]; mp = prtmas[ip];
   r1[2] = genypo[ig]; p1[2] = genymo[ig]; dp = gendmo[ig];
   r1[3] = genzpo[ig]; p1[3] = genzmo[ig]; r1[0] = 0.0;
   if(dp > 0.0)
   {
      pm = sqrt(p1[1]*p1[1] + p1[2]*p1[2] + p1[3]*p1[3]); if(pm > 0.0)
      {
         dx = p1[1]/pm; dy = p1[2]/pm; dz = p1[3]/pm;
         pm = uti_rnn(pm,dp); if(pm < 0.0) pm = -pm;
         p1[1] = dx*pm; p1[2] = dy*pm; p1[3] = dz*pm;
      }
   }
   p1[0] = sqrt(mp*mp + p1[1]*p1[1] + p1[2]*p1[2] + p1[3]*p1[3]);
   trknum++; it = trknum - 1;
   trktpo[it][0] = r1[0]; trkemo[it][0] = p1[0]; trkpnt[it] = 1;
   trkxpo[it][0] = r1[1]; trkxmo[it][0] = p1[1]; trkprt[it] = ip;
   trkypo[it][0] = r1[2]; trkymo[it][0] = p1[2];
   trkzpo[it][0] = r1[3]; trkzmo[it][0] = p1[3];
   if(genflg[ig] != 0)
   {
      id = abc_dec(it,genflg[ig]); it++;
   }

/* Swim tracks */
   while(it < trknum)
   {
   /* Initialize particle, step */
      ip = trkprt[it];
      is = 1;
      idt = 0;
      while(is < maxpnt)
      {
      /* Fill initial point */
         io = 0;
         id = 0;
         r1[0] = trktpo[it][is-1]; p1[0] = trkemo[it][is-1];
         r1[1] = trkxpo[it][is-1]; p1[1] = trkxmo[it][is-1];
         r1[2] = trkypo[it][is-1]; p1[2] = trkymo[it][is-1];
         r1[3] = trkzpo[it][is-1]; p1[3] = trkzmo[it][is-1];

      /* Determine initial object */
         if(is == 1)
         {
            i1 = -1; for(i=0;i<maxobj;i++) if(i == 0 || objnum[i] > 0)
            {
               if(uti_loc(i,r1[1],r1[2],r1[3]) < 1) i1 = i;
            }
         }

      /* Determine step direction,size */
         pm = sqrt(p1[1]*p1[1] + p1[2]*p1[2] + p1[3]*p1[3]);
         if(pm < minmom)
         {
            is += maxpnt;
            goto endstep;
         }
         dt = p1[0]/pm; dx = p1[1]/pm; dy = p1[2]/pm; dz = p1[3]/pm;
         rm = maxstp; if(prtchr[ip] != 0.0)
         {
            if(magtag && objbmg[i1] != 0.0)
            {
               rb = incmag*pm/(0.3*prtchr[ip]*objbmg[i1]);
               if(rb < 0) rb = -rb;
               if(rb >= minstp && rb < rm) rm = rb;
            }
            if(lostag && matrad[objmat[i1]] > 0.0)
            {
               rl = 0.01*incrad*matrad[objmat[i1]];
               if(rl >= minstp && rl < rm) rm = rl;
            }
         }

      /* Determine step point */
         r2[1] = r1[1] + rm*dx;
         r2[2] = r1[2] + rm*dy;
         r2[3] = r1[3] + rm*dz;
         /* Initialize object 2 */
         i2 = -1; for(i=0;i<maxobj;i++) if(i == 0 || objnum[i] > 0)
         {
            if(uti_loc(i,r2[1],r2[2],r2[3]) < 1) i2 = i;
         }
         /* Check if exit object 1 */
         ir = uti_int(i1,2,r1[1],r1[2],r1[3],r2[1],r2[2],r2[3],&x,&y,&z);
         if(ir == -1 || ir == +2)
         {
            io = -1; r2[1] = x; r2[2] = y; r2[3] = z;
            if(i2 > -1) for(i=0;i<maxobj;i++) if(i!=i1 && (i==0||objnum[i]>0))
            {
               if(uti_loc(i,r2[1],r2[2],r2[3]) < 1) i2 = i;
            }
            if(trctag == 10)
            fprintf(stderr,"crs2: %d %d %d %d %lf %lf %lf \n",
               is,i1,i2,ir,x,y,z);
         }
         /* Check if enter new object */
         for(i=1;i<maxobj;i++) if(i!=i1 && (i==0||objnum[i]>0))
         {
            ir = uti_int(i,1,r1[1],r1[2],r1[3],r2[1],r2[2],r2[3],&x,&y,&z);
            if(ir > 0)
            {
               io = +1; i2 = i; r2[1] = x; r2[2] = y; r2[3] = z;
               if(trctag == 10)
               fprintf(stderr,"crs1: %d %d %d %d %lf %lf %lf \n",
                  is,i1,i2,ir,x,y,z);
            }
         }
         /* Check for minimum step */
         rm = sqrt((r2[1]-r1[1])*(r2[1]-r1[1])+(r2[2]-r1[2])*(r2[2]-r1[2])
                  +(r2[3]-r1[3])*(r2[3]-r1[3]));
         if(rm < minstp)
         {
            if(trctag == 10) fprintf(stderr,"abc_gen: minimum step taken \n");
            rm = minstp; i2 = i1; io = -2;
            r2[1] = r1[1] + rm*dx;
            r2[2] = r1[2] + rm*dy;
            r2[3] = r1[3] + rm*dz;
         }
         /* Calculate step time,momentum */
         r2[0] = r1[0] + rm*dt;
         p2[0] = p1[0]; p2[1] = p1[1]; p2[2] = p1[2]; p2[3] = p1[3];

      /* Fill step point */
         trktpo[it][is] = r2[0]; trkemo[it][is] = p2[0];
         trkxpo[it][is] = r2[1]; trkxmo[it][is] = p2[1];
         trkypo[it][is] = r2[2]; trkymo[it][is] = p2[2];
         trkzpo[it][is] = r2[3]; trkzmo[it][is] = p2[3];
         trkpnt[it]++;
         if(trctag == 10)
         {
            fprintf(stderr,"gen: %d %d %d   %d %d   %lf %lf   %lf %lf \n",
               it,ip,is,i1,i2,r1[1],r1[2],r2[1],r2[2]);
         }

      /* Apply physical processes */
         if(magtag) if(prtchr[ip] != 0.0) abc_mag(it,i1,i2);    /* mag field */
         if(mlttag) if(prtchr[ip] != 0.0) abc_mlt(it,i1,i2);    /* mlt scatt */
         if(lostag) if(prtchr[ip] != 0.0) abc_los(it,i1,i2);  /* energy loss */
         if(dectag) if(prtmdn[ip] > 0) id = abc_dec(it,0);          /* decay */

      /* Update current object */
         if(io == 0)
         {
            i3 = -1; for(i=0;i<maxobj;i++) if(i == 0 || objnum[i] > 0)
            {
               if(uti_loc(i,r2[1],r2[2],r2[3]) < 1) i3 = i;
            }
            if(i3 > i2)
            {
               if(trctag == 1)
               fprintf(stderr,"gen: enter %d from %d \n",i3,i2);
               i2 = i3;
            }
            if(i3 < i2)
            {
               if(trctag == 1)
               fprintf(stderr,"gen: exit %d to %d \n",i2,i3);
               i2 = i3;
            }
         }

      /* Update current position,momentum */
         r2[0] = trktpo[it][is]; p2[0] = trkemo[it][is];
         r2[1] = trkxpo[it][is]; p2[1] = trkxmo[it][is];
         r2[2] = trkypo[it][is]; p2[2] = trkymo[it][is];
         r2[3] = trkzpo[it][is]; p2[3] = trkzmo[it][is];
         pm = sqrt(p2[1]*p2[1] + p2[2]*p2[2] + p2[3]*p2[3]);

      /* Check limits */
         if(i2 < 0) is += maxpnt;
         if(pm < minmom) is += maxpnt;
         if(id == 1) is += maxpnt;

      /* Handle detector output */
         if(idt == 0 && objdet[i1] > 0 && is ==  1)
         {
            idt = i1;
            r1d[0] = r1[0]; r1d[1] = r1[1]; r1d[2] = r1[2]; r1d[3] = r1[3];
            p1d[0] = p1[0]; p1d[1] = p1[1]; p1d[2] = p1[2]; p1d[3] = p1[3];
         }
         if(idt == 0 && objdet[i2] > 0 && i2 != i1)
         {
            idt = i2;
            r1d[0] = r2[0]; r1d[1] = r2[1]; r1d[2] = r2[2]; r1d[3] = r2[3];
            p1d[0] = p2[0]; p1d[1] = p2[1]; p1d[2] = p2[2]; p1d[3] = p2[3];
         }
         if(idt != 0 && objdet[i1] > 0 && (i2!=i1 || is>=(maxpnt-1)))
         {
            r2d[0] = r2[0]; r2d[1] = r2[1]; r2d[2] = r2[2]; r2d[3] = r2[3];
            p2d[0] = p2[0]; p2d[1] = p2[1]; p2d[2] = p2[2]; p2d[3] = p2[3];
            if(dattag == 1)
            {
               fprintf(stdout,"%d %d %d ",idt,ig,ip);
               for(i=0;i<4;i++) fprintf(stdout,"%f ",r1d[i]);
               for(i=0;i<4;i++) fprintf(stdout,"%f ",p1d[i]);
               for(i=0;i<4;i++) fprintf(stdout,"%f ",r2d[i]);
               for(i=0;i<4;i++) fprintf(stdout,"%f ",p2d[i]);
               fprintf(stdout,"\n");
            }
            if(dattag == 2)
            {
               fprintf(stdout,"  %02d %02d %02d  ",idt,ig,ip);
               for(i=0;i<4;i++) fprintf(stdout,"%7.3f ",r1d[i]);
               fprintf(stdout,"  ");
               for(i=0;i<4;i++) fprintf(stdout,"%7.3f ",p1d[i]);
               fprintf(stdout,"\n"); fprintf(stdout,"            ");
               for(i=0;i<4;i++) fprintf(stdout,"%7.3f ",r2d[i]);
               fprintf(stdout,"  ");
               for(i=0;i<4;i++) fprintf(stdout,"%7.3f ",p2d[i]);
               fprintf(stdout,"\n");
            }
            idt = 0;
            if(hitnum < maxhit)
            {
               hitnum++; ih = hitnum - 1; hitprt[ih] = ip;
               hitxpo[ih] = r1d[1]; hitypo[ih] = r1d[2]; hitzpo[ih] = r1d[3];
               hitnum++; ih = hitnum - 1; hitprt[ih] = ip;
               hitxpo[ih] = r2d[1]; hitypo[ih] = r2d[2]; hitzpo[ih] = r2d[3];
            }
            else fprintf(stderr,"  abc_gen: maximum hits reached \n");
         }

      /* Increment object,step */
         endstep:
         i1 = i2;
         is++;
      }
      it++;
   }
}