*	This is a simple test file which initialize muon cross-sections
*	and energy losses (initialize_music) and propagate muon down to
*	a predefined depth (muon_transport) in standard rock and sea water. 
*       It uses the following files:
*	- music-*-2004.f - file with all subroutines including 
*		initialize_music, muon_transport, music etc.;
*	- music-cross-sections-*.dat - formatted files with tables with 
*		integral cross-sections;
*	- music-eloss-*.dat - formatted files with continuous muon energy
*		losses (v<v_cut=10^(-3));
*	- music-double-diff-*.dat - formatted files with inelastic
*		scattering cross-section to sample angular deviation.
*
*	Initial parameters:
*	xi,yi,zi - initial muon coordinates in cm;
*	cxi,cyi,czi - initial directional cosines;
*	emui - muon energy in GeV;
*	depth* - distance to the next output from the code (distance to propagate
*	muon) 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;
*       
*	Program returns the values of variables after the muon was 
*       propagated down to 'depth'.
*	Output parameters:
*	xf,yf,zf - final muon coordinates in cm;
*	cxf,cyf,czf - final directional cosines;
*	emuf - muon energy in GeV;
*	depthf - distance to which the muon was propagated in cm 
*       (should be similar to 'depth0')
*	ttimef - time in nsec;
*
*!       An input parameter 'ic' is passed to the subroutine
*!       initialize_music_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 muons 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 muon energy is 10^4 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_music_sr(ic)  ! initialisation of cross-sections for standard rock
	call initialize_music_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 withe respect to z-axis
	emu0=1.d4  ! 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 muons in standard rock in cm
	depth_water=50000.  ! distance to transport muons 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.

*       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

*       Transport in standard rock

	call muon_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)

*       Transport in sea water

	call muon_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)

*       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

	if(emuf.gt.0.106) 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
	end if
c         print *,'theta0,phi0  ',theta0,phi0
c         print *,'theta,phi,  ',theta,phi
c         print *,'x0,y0,z0  ',x0,y0,z0
c         print *,'x,y,z  ',x,y,z
c         print *,'emu0  ',emu0
c         print *,'emu  ',emu		!if emu<0.106 muon was stopped
      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
      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

      stop
      end