*	This is a simple test file which initialize tau cross-sections
*	and energy losses (initialize_tausic) and propagate tau to
*	a predefined distance (tau_transport) in standard rock and sea water. 
*       It uses the following files:
*	- tausic-*-2004.f - file with all subroutines including 
*		initialize_tausic, tau_transport, tausic etc.;
*	- tausic-cross-sections-*.dat - formatted files with tables with 
*		integral cross-sections;
*	- tausic-eloss-*.dat - formatted files with continuous tau energy
*		losses (v<v_cut=10^(-3));
*	- tausic-double-diff-*.dat - formatted files with inelastic
*		scattering cross-section to sample angular deviation.
*
*	Initial parameters:
*	xi,yi,zi - initial tau coordinates in cm;
*	cxi,cyi,czi - initial directional cosines;
*	emui - tau energy in GeV;
*	depth* - distance to the next output from the code (distance to propagate
*	tau) in cm;
*	tmui - time in nsec;
*       rho* - density of medium; default: rho_rock = 2.65 for rock, 
*       rho_water = 1.03 for sea water;
*       ms_flag - flag to switch on/off multiple scattering: 1 - on, 0 - off;
*       d_flag - flag to switch on/off angular deflection due to other processes: 
*       1 - on, 0 - off;
*       itflag - flag to identify whether the (remaining) life time of tau 
*                (in tau coordinate system) was calculated or not: 
*                itflag = 0 - not calculated (start of the tau track);
*                itflag = 1 - calculated (transfered from previous call to transport);
*	tauti - initial (remaining) life time of tau (should be 0 if itflag = 0).
*       
*	Program returns the values of variables after the tau was 
*       propagated down to 'depth'.
*	Output parameters:
*	xf,yf,zf - final tau coordinates in cm;
*	cxf,cyf,czf - final directional cosines;
*	emuf - tau energy in GeV;
*	depthf - distance to which the tau was propagated in cm; 
*       (should be similar to 'depth0')
*	ttimef - time in nsec;
*       idec - flag specifying whether a tau decayed (1) or not (0);
*	tautf - final (remaining) life time of tau - output from propagation.
*
*!       An input parameter 'ic' is passed to the subroutine
*!       initialize_tausic_sr(ic)
*!       If ic = 1 then the inelastic cross-section calculated by
*!       Bugaev&Shlepin is used
*!       If ic = 2 then the inelastic cross-section calculated using
*!       structure functions is used
*
*	Note: the subroutines use standard CERN library routines and 
*	random number generators RANLUX. The code uses also
*	the generator of a pair of normally distributed numbers with a 
*	predefined correlation coefficient CORSET and CORGEN
*	from CERN library. The sequence of random numbers should be
*	initialized at the beginning of each run using RLUXGO and
*	RMARIN (see CERN library documentation)
*
*	As an example, this program propagates 1000 taus to a distance of
*	1000 m in standard rock (note: conversion from 
*	cm to m w.e. is done inside the code) and then to a distance of
*       500 m in sea water.
*       Initial tau energy is 10^9 GeV
*
	implicit real*8 (a-h,o-z)
	parameter (pi=3.141592654d0)
	real*4 yfl
	integer*4 ms_flag,d_flag

*       Initialisation

        ic = 1 ! inelastic cross-section from Bugaev&Shlepin
*        ic = 2 ! inelastic cross-section with structure functions
	call initialize_tausic_sr(ic)  ! initialisation of cross-sections for standard rock
	call initialize_tausic_sw(ic)  ! initialisation of cross-sections for sea water
	iranlux=1  ! initial seed for ranlux
	iranlux1=1 ! initial seed for ranmar
	call rluxgo(3,iranlux,0,0)  ! initialisation of ranlux
	call rmarin(iranlux1,0,0)   ! initialisation of ranmar
	nmumax=1000 ! number of muons to transport
	theta0=0.   ! initial zenith angle
	phi0=0.     ! initial azimuthal angle
	x0=0.       ! initial x
	y0=0.       ! initial y
	z0=0.       ! initial z - theta is calculated with respect to z-axis
	emu0=1.d9   ! initial muon energy
	rho_rock = 2.65  ! density for standard rock
	rho_water = 1.03 ! density for sea water
	depth_rock=100000.  ! distance to transport taus in standard rock in cm
	depth_water=50000.  ! distance to transport taus in sea water in cm
	ms_flag = 1         ! flag to switch on (1) and off (0) multiple scattering
	d_flag = 1          ! flag to switch on (1) and off (0) deflection due to stochastic processes
	ttime0=0.           ! initial time

	kmu=0      ! muon counter
	eemu=0.    ! this and next few lines initialise counters to calculate mean muon
	rmu=0.     ! parameters (useful for tests only)
	xmu=0.
	ymu=0.
	zmu=0.
	themu=0.
	phimu=0.
	timu = 0.
	dec = 0

*       Transport

	do i=1,nmumax          ! loop to transport muons
	emui=emu0   ! all initial muon energies and other parameters are defined above
	xi=x0
	yi=y0
	zi=z0
	theta=theta0
	phi=phi0
	cz0=dcos(theta0)               ! calculation of directional cosines
	cx0=dsin(theta0)*dcos(phi0)
	cy0=dsin(theta0)*dsin(phi0)
	cxi=cx0
	cyi=cy0
	czi=cz0
	emui=emu0
	tmui=ttime0

	idec = 0    ! output parameter for tau decay; idec = 0 (no decay), idec = 1 (tau decayed)
	itflag = 0  ! the (remaining) life time of tau (in tau coordinate system) was not calculated yet
	tauti = 0.  ! initial (remaining) life time of tau (should be 0 if itflag = 0)
	tautf = 0.  ! final (remaining) life time of tau - output from propagation

*       Transport in standard rock

	call tau_transport_sr(xi,yi,zi,cxi,cyi,czi,emui,depth_rock,
     *  tmui,rho_rock,ms_flag,d_flag,xf1,yf1,zf1,cxf1,cyf1,czf1,Emuf1,
     *  depthf1,tmuf1,idec1,itflag,tauti,tautf1)

c	call tau_transport_sr(xi,yi,zi,cxi,cyi,czi,emui,depth_rock,
c     *  tmui,rho_rock,ms_flag,d_flag,xf,yf,zf,cxf,cyf,czf,Emuf,
c     *  depthf,tmuf,idec,itflag,tauti,tautf)

*       Transport in sea water

c	call tau_transport_sw(xi,yi,zi,cxi,cyi,czi,emui,depth_water,
c     *  tmui,rho_water,ms_flag,d_flag,xf,yf,zf,cxf,cyf,czf,Emuf,
c     *  depthf,tmuf,idec,itflag,tauti,tautf)

	xf = xf1
	yf = yf1
	zf = zf1
	cxf = cxf1
	cyf = cyf1
	czf = czf1
	Emuf = Emuf1
	tmuf = tmuf1
	tautf = tautf1
	idec = idec1

	if(idec1.eq.0) then

	   tauti = tautf1
	   itflag = 1
	   idec = 0
	call tau_transport_sw(xf1,yf1,zf1,cxf1,cyf1,czf1,emuf1,
     *  depth_water,tmuf1,rho_water,ms_flag,d_flag,xf,yf,zf,
     *  cxf,cyf,czf,Emuf,depthf,tmuf,idec,itflag,tauti,tautf)

	end if

*       Calculation of theta and phi - useful for tests only

	theta=dacos(czf)
	IF(cxf.NE.0.) THEN
	phi=datan(cyf/cxf)
	IF(cxf.LT.0.) phi=pi+phi
	END IF
	IF(cxf.EQ.0.) THEN
	phi=pi/2.
	IF(cyf.LT.0.) phi=pi*3./2.
	END IF
	IF(phi.GT.pi*2.) phi=-pi*2.+phi
	IF(phi.LT.0.) phi=phi+pi*2.
	theta = theta * 180. / pi
	phi = phi * 180. / pi

*       Summing up muon parameters for calculation of average values - useful for tests only

c	if(emuf.gt.1.778.and.idec.eq.0) then
	if(emuf.gt.1.778) then
	kmu=kmu+1
	eemu=eemu+emuf
	xmu=xmu+xf
	ymu=ymu+yf
	zmu=zmu+zf
	themu=themu+theta
	phimu=phimu+phi
	rmu=sqrt(xf*xf+yf*yf)+rmu
	timu = timu + tmuf
	dec = dec + idec
	end if
c	if(emuf.le.1.778) then
c         print *,'theta0,phi0  ',theta0,phi0
c         print *,'theta,phi,  ',theta,phi
c         print *,'x0,y0,z0  ',xi,yi,zi
c         print *,'x,y,z  ',xf,yf,zf
c         print *,'emu0  ',emu0
c         print *,'emu  ',emuf		!if emu<1.778 tau was stopped
c         print *,'tmuf  ',tmuf
c	 print *,'idec ',idec
c	end if
	end do

*       Calculation of average values for muon parameters after transport

      if(kmu.gt.0) then
      spmu=1.*kmu/nmumax
      eemu=eemu/kmu
      xmu=xmu/kmu
      ymu=ymu/kmu
      zmu=zmu/kmu
      themu=themu/kmu
      phimu=phimu/kmu
      rmu=rmu/kmu
	timu = timu / kmu
	dec = dec / kmu
      end if
      print *,'survival probability',spmu
      print *,'mean energy in GeV',eemu
      print *,'x,y,z in cm',xmu,ymu,zmu
      print *,'lateral displacement in cm',rmu
      print *,'theta in degrees',themu
      print *,'phi in degrees',phimu
	print *,'mean time',timu
	print *,'mean decay flag',dec

      stop
      end