! program to test injection line field maps. ! In the lattice field, the fields are scaled by 1/2 so that the on axis trajectory stays within the boundary of the grids ! In the program, the field is restored to full scale and the trajectory is initialized so that it can stay with the grid map ! Computes trajectory and field along the trajectory and writes to fort.14 ! Computes field on axis and writes to file fort.13 ! ! Reads lattice file g-2/files/bmad_grid_fringe_inf.bmad ! The fringe and inflector grid files are in g-2/files/backleg/ ! Latest version has FRINGE and INFLECTOR superimposed in a single grid file. ! The goal is to separate so that the inflector orientation and aperture can be manipulated independently ! program test_em_field_calc use bmad use muon_mod use random_mod use materials_mod use parameters_bmad implicit none type (lat_struct) lat type (coord_struct) co, orb_start, orb_end, orb_init, orb type (coord_struct), allocatable :: orbit(:) type (em_field_struct) field, field_mpole type(ele_struct) ele type (ele_struct), target :: ele_p type (ele_struct), pointer :: lord type (muon_struct), allocatable :: muons(:), muons_ele(:,:) integer lun integer lun1 integer nargs, iargc, ios integer ix integer track_state integer i integer n integer j integer nmuons integer, allocatable :: nalive(:) integer status/0/ integer datetime_values(8) integer system character*120 line, lat_file, lat_file_name character*120 string character*3 cmax_loops character*16 date, time, zone character*10 snumber logical err_flag, local_ref_frame/.false./ logical grid logical s_to_s logical create_a_distribution/.false./ real(rp) y, s_pos, x real(rp) s_next real(rp) vec(6) real(rp) sigma_x,epsx,sigma_xp real(rp) betax, alphax, betay, alphay call date_and_time(date,time,zone,datetime_values) print '(a10,a10,a1,a10)',' date = ',date,' time = ', time directory = trim(date)//'_'//trim(time(1:6)) status= system('mkdir ' // directory) if(status == 0)print '(a20,a)',' create directory = ', trim(directory) us_scatter= .false. ds_scatter = .false. call initializeMaterials() inflector_width = 0 nmuons = 100 nargs = command_argument_count() if (nargs >= 1) call get_command_argument(1, lat_file) if (nargs == 2) then call get_command_argument(2, snumber) read(snumber,*)nmuons endif if(nargs == 0)then lat_file = 'bmad.' print '(a,$)',' Lattice file name? (Default = bmad.) ' read(5,'(a)') line call string_trim(line, line, ix) lat_file = line if (ix==0) lat_file = 'bmad.' endif print '(a,a,a,1x,i10)', ' lat_file = ',trim( lat_file), ' nmuons = ', nmuons string = 'cp'//' '//trim(lat_file)//' '//trim(directory)//'/.' print *,string call execute_command_line(string) call bmad_parser(lat_file, lat, make_mats6 = .false.) call reallocate_coord(orbit, lat%n_ele_max) orbit(0)%vec=0 orbit(0)%vec(1) = -0.047671-0.00005 orbit(0)%vec(2) = .022134-0.00005 !0.021 orbit(0)%vec(3) = 0.0 !orbit(0)%vec(6) = -0.9 call init_coord(orb_start) orb_start = orbit(0) orb_init = orbit(0) do j=1, lat%n_ele_max ele_p = lat%ele(j) grid=.false. if((associated(ele_p%grid_field)))then grid = .true. lord => ele_p elseif(ele_p%field_calc == refer_to_lords$)then do i=1,ele_p%n_lord lord => pointer_to_lord(ele_p,i) if(associated(lord%grid_field))then grid = .true. exit endif end do endif if(grid)then print *,' grid = ', grid do i = 1, size(lord%grid_field) lord%grid_field(i)%field_scale = 1. end do endif end do !lat%ele(1)%grid_field(1)%field_scale = 1. !lat%ele(1)%grid_field(2)%field_scale = 1.012747 print *,' call track_all' call track_all(lat, orbit) do i=1,lat%n_ele_track print '(7es12.4)',lat%ele(i)%s,orbit(i)%vec write(15,'(7es12.4)') lat%ele(i)%s,orbit(i)%vec end do s_to_s=.true. if(s_to_s)then s_pos = 0 print '(a,6es12.4)' ,' start s to s tracking: orb_start%vec = ',orb_start%vec ix=1 do while(s_pos + 0.01 < lat%ele(lat%n_ele_track)%s) s_next = s_pos + 0.01 do i=1,lat%n_ele_track-1 if(lat%ele(i-1)%s 4.2)then !scan x-y plane for B orb_end%vec=0 orb_end%vec(3)=-0.028 do while( orb_end%vec(3) < 0.028) orb_end%vec(3) = orb_end%vec(3)+0.001 orb_end%vec(1) =-0.009 do while(orb_end%vec(1) < 0.009 .and. orb_end%vec(3) < 0.028) orb_end%vec(1) = orb_end%vec(1)+0.001 call em_field_calc(lat%ele(ix), lat%param, s_next-lat%ele(ix-1)%s, orb_end, local_ref_frame,field,err_flag=err_flag) write(16,'(10es12.4)')s_next, orb_end%vec, field%b write(6,'(10es12.4)')s_next, orb_end%vec, field%b end do end do endif write(14,'(10es12.4)')s_next, orb_end%vec, field%b s_pos = s_next end do endif print *,' track s to s complete' us_scatter = .false. ds_scatter = .false. sigma_x = 0.001 epsx=10.e-6 sigma_xp = epsx/sigma_x betax=50 alphax=-16 betay=10 alphay=-3 if(create_a_distribution)then !create a distribution allocate(muons(nmuons)) allocate(muons_ele(nmuons,1:lat%n_ele_track)) allocate(nalive(1:lat%n_ele_track)) nalive(:) = 0 call ran_gauss(muons(:)%gas) muons(:)%coord%vec(1) = muons(:)%gas*sqrt(epsx*betax) call ran_gauss(muons(:)%gas) muons(:)%coord%vec(2) = muons(:)%gas*sqrt(epsx/betax)*alphax call ran_gauss(muons(:)%gas) muons(:)%coord%vec(3) = muons(:)%gas * sqrt(epsx*betay) call ran_gauss(muons(:)%gas) muons(:)%coord%vec(4) = muons(:)%gas * sqrt(epsx*betay)*alphay do j=1,nmuons orbit(0)%vec = muons(j)%coord%vec + orb_init%vec !print *,' call track_all distribution' call track_all(lat, orbit) do i=1,lat%n_ele_track if(orbit(i)%state == alive$)nalive(i) = nalive(i) + 1 muons_ele(j,i)%coord = orbit(i) end do end do lun=lunget() open(unit=lun,file=trim(directory)//'/'//'total_muons.dat') lun1 = lunget() open(unit=lun1,file=trim(directory)//'/'//'phase_space.dat') do i=1,lat%n_ele_track write(lun,'(i10,1x,es12.4,1x,i10)')i, lat%ele(i)%s, nalive(i) write(lun1,'(/,a1,1x,i10,1x,es12.4,/)')'#',i, lat%ele(i)%s do j=1,nmuons if(muons_ele(j,i)%coord%state == alive$)write(lun1,'(7es12.4)')lat%ele(i)%s,muons_ele(j,i)%coord%vec end do end do close(unit=lun) close(unit=lun1) endif !create a distribution if(status == 0)print '(a20,a)',' create directory = ', trim(directory) end