C+++
C	PROGRAM		btraj
C
C	PURPOSE		To compute the trajectory of an electron in a magnetic
C			field B (read from a file)		
C
C	ALGORITHM	See formulae in ESRF red book	pag CIV-297
C
C	COMMENTS	Uses SHADOW's reference frame
C
C	CREATION DATE	10/91 M. Sanchez del Rio, C. Vettier
C
C--
	PROGRAM		btraj
	IMPLICIT REAL*8	(A-H,O-Z)
	PARAMETER	(nmax=5001)   !imposed by NPHOTON
	CHARACTER *80 	RSTRING,OUTFILE,INFILE
	DIMENSION 	BETAX(nmax), BETAZ(nmax),BETAY(nmax)
	DIMENSION	YX(nmax),YZ(nmax),CURV(nmax),YY(nmax)
	DIMENSION	BX(nmax),BZ(nmax),BH(10)
c
     	DATA	PI     	/  3.1415 92653 58979 32384 62643 D0 /
     	DATA	PIHALF 	/  1.5707 96326 79489 66192 31322 D0 /
     	DATA	TWOPI 	/  6.2831 85307 17958 64769 25287 D0 /
     	DATA	TODEG 	/ 57.2957 79513 08232 08767 98155 D0 /
     	DATA	TORAD	/  0.0174 53292 51994 32957 69237 D0 /
	DATA	TOCM	/  1.239 852	D-4		     /
	DATA	TOANGS 	/  1.239 852    D+4		     /
c
        c  = 2.998D8            !speed of light, m/s
        rm = 9.109D-31          !electron rest mass  kg
        e  = 1.602D-19          !electron charge, C
        h  = 6.626D-34          !Planck's constant   joules*sec
        c2 = e/(TWOPI*rm*c)     !unit(coul*sec/(kg*m))
c
	type *,'Wiggler trajectory calculation from: '
	type *,'  [0] An input file (2-column) with the magnetic field',
     $         '     of one period. 1st column: y[m]; 2nd column: B[T]'
	type *,'  [1] A file with magnetic field hatmonic data:',
     $         '    1st column: # (harmonic number); 2nd column: B[T]' 
	accept *,i_file
c
	infile = rstring ('Input file ? ')
	open (21,file=infile,status='old')
c
c	reads file with magnetic field Bz[T](y[m]) (vertical)
c
	if (i_file.eq.0) then
	 do i=1,nmax
	  read (21,*,end=33) YY(i),BZ(i)
	 end do
33 	 continue
	 npoints = i-1
	 ystep = yy(2)-yy(1)
	 per  = yy(npoints)-yy(1)
	 type *,npoints ,' points read from file'
	 type *,' Period of the ID is : ',per
	 nper = irint ('Enter number of periods : ')
c
c	read file with period [m], number of periods, harmonics [T]
c
	else if (i_file.eq.1) then
	 per = rnumber('Enter wiggler wavelength [m] : ')
	 nper = irint ('Enter number of periods : ')
c	get harmonics from file
	 do i=1,3001
	  read (21,*,end=34) n,bh(n)
	 end do
34 	 continue
	 nharm = i-1
	 type *,nharm,' harmonics read from file'
c	 type *,' Period of the ID is : ',per ,' [m]'
c	 type *,' Number of periods : ',nper
	 type *,'Enter number of points per period : '
	 accept *,npoints
	 ystep = per/(npoints-1)
	end if
	close (21)
	if (npoints*nper.gt.nmax) type *,'Too many points. Reenter file with
     $ less points [max=5001]'
c
c	calculates initial electron speed from its energy
c
	oldener = rnumber  ('Electron energy [GeV] ? ')
	ENER    = E*oldENER*1.0D9	!CHANGE UNITS (TO JOULE)
	GA0     = ENER/RM/C**2		!Gamma_0, in electron rest mass
     					!energy (mc^2) units,
     					!adimensional
	emc = e/(ga0*rm*c)
	beta02 = 1.0d0 - (rm*c**2/ener)**2
	beta0 = sqrt(beta02)
	type *,' beta = v/c =',beta0
	start_len=per*(nper-1)/2.d0		
c
c	calculates the integral of the magnetic field bz (proportional
c	to speeds)
c	
	if (i_file.eq.0) then
	 yint = 0.0d0
	 betax(1) = 0.0d0
	 do i=2,npoints
	   yint = yint + bz(i)
	   betax (i) = yint*ystep
	 end do 
	else if (i_file.eq.1) then
	 phase0=-pi
	 do i=1,npoints
		yy(i)=(i-1)*ystep-per/2.d0
		bz(i)=0.d0
		betax(i)=0.d0
		do n=1,nharm
			phase=twopi*(yy(i)/per)
			bz(i)=bz(i)+bh(n)*dcos(phase*n)
			betax(i)=betax(i)+bh(n)
	1			*(dsin(phase*n)-dsin(phase0*n))/n
		enddo
		betax(i)=betax(i)*per/twopi
	 end do 
	end if
	type *,' Integral of B = ',betax(npoints), '  [Tm]'
c
c	rescale to speeds v/c = 0.3/E[GeV] * integral (Bz [T] ds[m])
c
	do i=1,npoints
	 betax (i) = - (0.3/oldener) * betax(i) 
	 betay (i) = sqrt(beta0**2 - betax(i)**2)
	 curv(i) = -emc*bz(i)/beta0
	end do	
c
c	calculates positions as the integral of speeds
c
	if (i_file.eq.0) then
	 yint = 0.0d0
	 do i=1,npoints
	  yint = yint + betax(i)
	  yx(i) = yint * ystep 
	 end do 
	else if (i_file.eq.1) then
	 do i=1,npoints
		yx(i)=0.d0
		do n=1,nharm
			phase=twopi*(yy(i)/per)
			yx(i)=yx(i)-bh(n)*
	1			(dcos(phase*n)-dcos(phase0*n))/n/n
		enddo
		yx(i)=yx(i)*(-3.d-1/oldener)*(per/twopi)**2
	 end do 
	end if
c
c	creates parameters file
c
	outfile = rstring ('Output file with parameters ? ')
	open(22,file=outfile,status='unknown')
	  write (22,*) '########### Data from btraj ##############'
	  write (22,*) 'Period of the ID is : ',per ,' [m]'
	  write (22,*) 'Number of periods : ',nper
	  write (22,*) 'Number of points per periods : ',npoints
	  write (22,*) 'beta = v/c =',beta0
	  write (22,*) 'Electron energy [GeV] = ',oldener
	if (i_file.eq.0) then
	  write(22,*) 'Magnetic field profile from file:   ',infile
          write(22,*) 'Integral of B = ',betax(npoints), '[Tm]'
	else if (i_file.eq.1) then 
	  write(22,*) 'Magnetic field harmonics from file ',infile
	  write(22,*) 'Number of harmonics: ',nharm
	end if
	close(22)
c
c	creates trajectory file
c
	outfile = rstring ('Output file with the trajectory ? ')
	open(22,file=outfile,status='unknown')
	do j=1,nper
		nn=1
		if(j.gt.1) nn=2 	! to avoid overlap
		do i=nn,npoints
		write (22,*) YX(i),yy(i)+(j-1)*per-start_len,0.0,
	1		betax(i),betay(i),0.0d0,curv(i)
		end do
	 end do
	close(22)
c
c	write B(y) in a file (if desired)
c
	if (i_file.eq.1) i_write =irint('Do you want to write the magnetic
     $ field in a file ? [0/1] ')
	if (i_write.eq.1) then
	 outfile = rstring ('Name of the file for B[T](y[m]) ? ')
	 open(22,file=outfile,status='unknown')
	 do j=1,nper
		nn=1
		if(j.gt.1) nn=2 	! to avoid overlap
		do i=nn,npoints
		write (22,*)  	yy(i)+(j-1)*per-start_len,bz(i)
		end do
	 end do
	close(22)
	end if
c
	type *,'All done'
	end