* ********************************************************* *
* Kumac for making colour plots out of FLUKA Usrbin outputs *
* ********************************************************* *
*                                                           *
* Initially written by:                                     *
*          Alfredo.Ferrari@cern.ch                          *
*          Paola.Sala@cern.ch                               *
*                                                           *
* Modified to plot also regions, materials, lattices and    *
* magnetic fields                                           *
* Vasilis.Vlachoudis@cern.ch 2004                           *
* ********************************************************* *
 opt *
 set *
 opt utit
 opt nbox

 close 1
 h/del *
 v/del *

* To remove the outer box
*opt nbox
 opt nstat
* Title font
 set *fon -60
* Axis label font
 set lfon -10
* Axis title size
 set asiz .5
* Simil bold
*set pass 5
*set cshi .01
* Distance of Y-axis numbers to the axis
 set xval .2
* Distance of X-axis numbers to the axis
 set yval .2
* Distance of X-axis label to the axis
 set ylab .6

* Global Variables
 FILE='pawlevh'
 FILEPG='plotgeom.his'

* Ask user
 read FILE 'Input file'
 psfile=[FILE]//'.ps'
 read psfile 'Postscript output file'
 shell rm [psfile]
 type=-112
 read type 'Metafile type (PS-PORT=-111, PS-LAND=-112, EPS=-113)'
 f/file 60 [psfile]
 meta 60 [type]
 set ygti .6
 set gsiz .45

*title
 read title 'Title'
 title [title]
 v/read xlim,ylim [FILE].lim
 v/read npoints [FILE].npo
 v/read x,y [FILE].poi

 message 'Input X/Y-axes stretching factor'
 factx=1.0
 read factx 'Input factor X'

 facty=1.0
 read facty 'Input factor Y'

*
 xmin=xlim(1)
 xmax=xlim(2)
 ymin=ylim(1)
 ymax=ylim(2)
 nc=30
 ncc=[nc]-7
 v/crea sizes(2) r
 xsiz=(xlim(2)-xlim(1))*[factx]
 ysiz=(ylim(2)-ylim(1))*[facty]
 v/inp sizes [xsiz] [ysiz]
 sigma sm=vmax(sizes)
 
* 19 cm width or height (the longest)
 fact=19./sm(1)
 xsiz=[xsiz]*[fact]
 ysiz=[ysiz]*[fact]
 yasz=0.3/(.76*[ysiz])

 if [xsiz]<11 then
   xsiz=0.76*[xsiz]+1.54+1.10
   size [xsiz] [ysiz]
   set xmgr 1.54
   set xmgl 1.10
   set ymgu .12*[ysiz]
   set ymgl .12*[ysiz]
 elseif [ysiz]<11 then
   ysiz=0.76*[ysiz]+1.32+1.32
   size [xsiz] [ysiz]
   set xmgr .14*[xsiz]
   set xmgl .10*[xsiz]
   set ymgu 1.32
   set ymgl 1.32
 else
   size [xsiz] [ysiz]
   set xmgr .14*[xsiz]
   set xmgl .10*[xsiz]
   set ymgu .12*[ysiz]
   set ymgl .12*[ysiz]
 endif

 igset mtyp 20
 igset lwid 2
 igset mscf .2
 igset fais 1
 set   asiz .4
 set   tsiz .5
 
* Select what to plot
 message 'Plot type:'
 message '	D=Data*'
 message '	P=R-Z-Phi'
 message '	G=Geometry only'
 message '	R=Regions'
 message '	M=Materials'
 message '	F=Magnetic Field'
 message '	L=Lattice'
 read plt
 plt=$upper([plt])

 if [plt]='G' then
   * Geometry Only
   message 'Geometry Plot'
   null [xmin] [xmax] [ymin] [ymax]
   exe plotgeom
   exe center

 elseif [plt]='R' then
   * Regions
   message 'Regions Plot'
   exe palette
   exe plothist 1 [FILEPG]
   exe plotgeom
   exe center

 elseif [plt]='M' then
   * Materials
   message 'Materials Plot'
   exe palette
   exe plothist 2 [FILEPG]
   exe plotgeom
   exe center

 elseif [plt]='L' then
   * Lattices
   message 'Lattices Plot'
   exe palette
   exe plothist 3 [FILEPG]
   exe plotgeom
   exe center

 elseif [plt]='F' then
   * Mangetic Field
   message 'Mangetic Field Plot'
   exe palette2
   exe plotfield [FILEPG]
   exe plotgeom
   exe legende
   exe center

 elseif [plt]='P' then
   * R-Z-Phi Plot
   message 'R-Z-Phi Plot'

   opt logz
   exe palette

   h/file 1 [FILE].his
   hrin 1

   nxh=$HINFO(1,'XBINS')
   nyh=$HINFO(1,'YBINS')
   xmih=$HINFO(1,'XMIN')
   xmah=$HINFO(1,'XMAX')
   ymih=$HINFO(1,'YMIN')
   ymah=$HINFO(1,'YMAX')
   xmin=-[ymah]
   xmax=[ymah]
   ymin=-[ymah]
   ymax=[ymah]
   read xmin
   read xmax
   read ymin
   read ymax

   null [xmin] [xmax] [ymin] [ymax]

   minval=xlim(5)
   maxval=ylim(5)
   message 'Min and Max value in matrix' [minval], [maxval]
   read minval
   read maxval

   *message enter n. of divisions per decade
   v/crea val(5) r
   v/crea nclim(3) i
   divn=3
   read divn 'Enter n. of divisions per decade'
   v/inp val [minval] [maxval] [divn] 0. 0.
   v/inp nclim 0 0 [ncc]
   call $FLUPRO/flutil/pawlevhlp.for(1.)
   minval=val(1)
   maxval=val(2)
   minsc=val(4)
   maxsc=val(5)
   nca=nclim(1)
   ncb=nclim(2)
   min 1 [minval]
   max 1 [maxval]
   
   ncuse=[ncb]-[nca]+1
   * call hgive(1,ct,nxh,xmih,xmah,nyh,ymih,ymah,nw,loc)
   * dp= ([xmah]-[xmih] )/[nxh]
   dp=360./[nxh]
   dr=([ymah]-[ymih] )/[nyh]
   v/crea ichan(3) i
   v/crea gival(4) r
   call $FLUPRO/flutil/fgivev.f(0)
   xmm=[xmax]+0.5*[dr]
   do ir=1,[nyh]
     rmi=( [ir] - 1 ) * [dr] + [ymih]
     rma = [rmi] + [dr]
     if [rma]<[xmm] then
       do ip = 1, [nxh]
	 pmi=( [ip] - 1 ) * [dp] -180.
	 pma = [pmi] + [dp]
	 igset fais 1
	 * con=$call('hij(1,[ip],[ir])')
	 v/inp ichan 1 [ip] [ir]
	 call $FLUPRO/flutil/fgivev.f(1)
	 icc=icval(1)
	 igset faci [icc]
	 arc 0. 0. [rmi] [rma] [pmi] [pma]
       enddo
     endif
   enddo
   
   * **************************************************************
   * colorbar
   dx=([xmax]-[xmin])/16
   xaxis=[xmax]+1.505*[dx]
   x2=[xaxis]-.005*[dx]
   x1=[xaxis]-[dx]
* igset laof ([ymax]-[ymin])/40.
* igset lasi ([ymax]-[ymin])/40.
   igset laof 2.8*[yasz]*([ymax]-[ymin])*[facty]
   igset lasi 0.85*[yasz]*([ymax]-[ymin])
   igset tmsi 0.7*[yasz]*([ymax]-[ymin])*[facty]
   axis [xaxis] [xaxis] [ymin] [ymax] [minsc] [maxsc] 211 'G=NHDS'
   dy=([ymax]-[ymin])/[ncuse]
   y2=[ymin]
   
   do ic=[nca],[ncb],1
     igset clip 0
     icc=7+[ic]
     y1=[y2]
     y2=[y1]+[dy]
     igset faci [icc]
     box [x1] [x2] [y1] [y2]
   enddo

   exe plotgeom
*
 else
   * Data
   opt logz
   exe palette
   
   h/file 1 [FILE].his
   minval=xlim(5)
   maxval=ylim(5)
   message 'Min and Max value in matrix' [minval], [maxval]
   read minval
   read maxval
   
   *message enter n. of divisions per decade
   v/crea val(5) r
   v/crea nclim(3) i
   divn=3
   read divn 'Enter n. of divisions per decade'
   v/inp val [minval] [maxval] [divn] 0. 0.
   v/inp nclim 0 0 [ncc]
   call $FLUPRO/flutil/pawlevhlp.for(1.)
   minval=val(1)
   maxval=val(2)
   minsc=val(4)
   maxsc=val(5)
   nca=nclim(1)
   ncb=nclim(2)
   min 1 [minval]
   max 1 [maxval]

   * uses nc-7 colors
   exe palette

   h/plo 1 colz([xmin]:[xmax],[ymin]:[ymax])

   exe plotgeom
 endif

 null [xmin] [xmax] [ymin] [ymax] SAB

 close 60
 close 1
 return

***************************************************************
* plot geometry
 macro plotgeom
   titxaxis='Z(cm)'
   message ' X axis label'
   read titxaxis
   tityaxis='R(cm)'
   message ' Y axis label'
   read tityaxis
   atit [titxaxis] [tityaxis] ! 333

   set plci 1
   sigma nv=nco(npoints)
   ilines=nv(1)
   ip=0
   ipa=1
   ipb=0
   do ip=1,[ilines]
     ipa=[ipb]+1
     np=npoints([ip])
     ipb=[ipa]-1+npoints([ip])
     if [np]=1 then
       pmarker [np] x([ipa]:[ipb]) y([ipa]:[ipb])
     elseif [np]>1 then
       pline  [np] x([ipa]:[ipb]) y([ipa]:[ipb])
     endif
   enddo
 return

***************************************************************
* plot a histogram with corrected colors
 macro plothist
   xmin=xlim(1)
   xmax=xlim(2)
   ymin=ylim(1)
   ymax=ylim(2)
   h/file 1 [2]
   hrin [1]
   call $FLUPRO/flutil/pawlevmat.for([1],22)
   MAX=$HINFO([1],'MAX')+1
   max [1] [MAX]
   h/plo [1] col([xmin]:[xmax],[ymin]:[ymax])
 return

***************************************************************
* plot magnetic field map
 macro plotfield
   h/file 1 [1]

   * Histograms
   HBX=4
   HBY=5
   HB=6

   * Plot limits
   XMIN=xlim(1)
   XMAX=xlim(2)
   YMIN=ylim(1)
   YMAX=ylim(2)

   SCALE=0.5
   ASCALE=0.1
   read SCALE 'Arrow length scaling factor'
   read ASCALE 'Arrow  size scaling factor'

   * Get field components in vectors
   hrin [HB]
   NX=$hinfo([HB],'XBINS')
   NY=$hinfo([HB],'YBINS')
   MINZ=$hinfo([HB],'MIN')
   MAXZ=$hinfo([HB],'MAX')
   if [MAXZ]=0 then
     null [xmin] [xmax] [ymin] [ymax]
     goto fin
   endif
   min [HB] 0.0
   h/plo [HB] colz([xmin]:[xmax],[ymin]:[ymax])

   N=[NX]*[NY]

   v/create bx([NX],[NY]) r
   v/create by([NX],[NY]) r
   h/get/cont [HBX] bx
   h/get/cont [HBY] by
   sigma b=sqrt(bx**2+by**2)

   XSTEP=([XMAX]-[XMIN])/[NX]
   YSTEP=([YMAX]-[YMIN])/[NY]

   X0=[XMIN]+[XSTEP]/2.0
   X1=[XMAX]-[XSTEP]/2.0
   sigma xp=array([NX],[X0]#[X1])

   Y0=[YMIN]+[YSTEP]/2.0
   Y1=[YMAX]-[YSTEP]/2.0
   sigma yp=array([NY],[Y0]#[Y1])

   if [NX]>50 then
     SX=INT([NX]/50)
   else
     SX=1
   endif

   if [NY]>50 then
     SY=INT([NY]/50)
   else
     SY=1
   endif

   do j=1,[NY],[SY]
     Y=[YMIN]+([j]-0.5)*[YSTEP]
     do i=1,[NX],[SX]
       X=[XMIN]+([i]-0.5)*[XSTEP]
       BX=bx([i],[j])
       BY=by([i],[j])
       B=b([i],[j])
       if [B]>0.00000001 then
	 X2=[X]+[BX]*[SCALE]
	 Y2=[Y]+[BY]*[SCALE]
	 SZ=[B]*[SCALE]*[ASCALE]
	 if [SZ]<0.001 then
	   SZ=0.001
	 endif

	 COL=1
	 if [BX]<0 then
	   COL=[COL]+2
	 endif
	 if [BY]<0 then
	   COL=[COL]+1
	 endif
	 set plci [COL]

	 arrow [X] [X2] [Y] [Y2] [SZ]

       endif
     enddo
   enddo
   v/del bx
   v/del by
   v/del b
fin:
 return

***************************************************************
* draw color bar
 macro colorbar
   do ic=1,[ncb]
     igset clip 0
     icc=7+[ic]
     y1=[y2]
     y2=[y1]+[dy]
     igset faci [icc]
     box [x1] [x2] [y1] [y2]
   enddo
 return

***************************************************************
* center cross
 macro center
   null 0 1 0 1 SAB
   set plci 1
   set lwid 1
   line 0.49 0.50 0.51 0.50
   line 0.50 0.49 0.50 0.51
 return

***************************************************************
* direction legende
 macro legende
   null 0 1 0 1 SAB

   XC=0.96
   YC=0.96
   R=0.03
   RA=0.02
   AZ=0.2

   set ltyp 21
   set plci 1
   set lwid 1
   set bord 1
   set fais 1
   set faci 0
   ellipse [XC] [YC] [R] [R]
   X1=[XC]-[R]
   X2=[XC]+[R]
   Y1=[YC]-[R]
   Y2=[YC]+[R]
   line [X1] [YC] [X2] [YC]
   line [XC] [Y1] [XC] [Y2]

   set ltyp 0
   set lwid 7

				| NE arrow
   set plci 1
   X=[XC]+[RA]
   Y=[YC]+[RA]
   arrow [XC] [X] [YC] [Y] [AZ]

				| SE arrow
   set plci 2
   X=[XC]+[RA]
   Y=[YC]-[RA]
   arrow [XC] [X] [YC] [Y] [AZ]

				| SW arrow
   set plci 4
   X=[XC]-[RA]
   Y=[YC]-[RA]
   arrow [XC] [X] [YC] [Y] [AZ]

				| NW arrow
   set plci 3
   X=[XC]-[RA]
   Y=[YC]+[RA]
   arrow [XC] [X] [YC] [Y] [AZ]
 return

***************************************************************
* color palette
 macro palette
   igset ncol 30

* Fluka colors using white as background
   color_table   30   0.0   0.0   0.0
   color_table   29   0.4   0.0   0.0
   color_table   28   0.6   0.0   0.0
   color_table   27   0.8   0.0   0.0
   color_table   26   1.0   0.0   0.0
   color_table   25   1.0   0.5   0.0
   color_table   24   1.0   0.8   0.0
   color_table   23   1.0   1.0   0.0
   color_table   22   0.8   1.0   0.0
   color_table   21   0.5   1.0   0.0
   color_table   20   0.0   0.9   0.2
   color_table   19   0.0   0.7   0.5
   color_table   18   0.0   0.8   1.0
   color_table   17   0.0   0.6   1.0
   color_table   16   0.0   0.0   1.0
   color_table   15   0.0   0.0   0.8
   color_table   14   0.5   0.0   0.8
   color_table   13   0.7   0.0   1.0
   color_table   12   0.9   0.0   1.0
   color_table   11   1.0   0.4   1.0
   color_table   10   0.9   0.6   0.9
   color_table    9   0.8   0.8   0.8
   color_table    8   1.0   1.0   1.0
 return

***************************************************************
* color palette2
 macro palette2
   igset ncol 30

* Fluka colors where COLOR = (1+COLOR)/2
   color_table   30   0.50  0.50  0.50
   color_table   29   0.70  0.50  0.50
   color_table   28   0.80  0.50  0.50
   color_table   27   0.90  0.50  0.50
   color_table   26   1.00  0.50  0.50
   color_table   25   1.00  0.75  0.50
   color_table   24   1.00  0.90  0.50
   color_table   23   1.00  1.00  0.50
   color_table   22   0.90  1.00  0.50
   color_table   21   0.75  1.00  0.50
   color_table   20   0.50  0.95  0.60
   color_table   19   0.50  0.85  0.75
   color_table   18   0.50  0.90  1.00
   color_table   17   0.50  0.80  1.00
   color_table   16   0.50  0.50  1.00
   color_table   15   0.50  0.50  0.90
   color_table   14   0.75  0.50  0.90
   color_table   13   0.85  0.50  1.00
   color_table   12   0.95  0.50  1.00
   color_table   11   1.00  0.70  1.00
   color_table   10   0.95  0.80  0.95
   color_table    9   0.90  0.90  0.90
   color_table    8   1.00  1.00  1.00
 return