*********************************************************************** SUBROUTINE GUFLDO(XV,HV) * XV: 3 dim Position vector (input) * HV: 3 dim Field vector (output) * * returns magnetic field at position XV * * last modified on 11/2/98 by R.Holzmann *********************************************************************** implicit none #include "user.inc" #include "fields.inc" REAL X_B, Y_B, Z_B, XV(3), HV(3) REAL FAC, PI, PHI, PHI_GAP, PHI_GAP_ABS, HX, HY, HZ REAL*4 CPHI, SPHI, PHI0, WEIGHT_Z, WEIGHT_RHO, WEIGHT_PHI, RHO REAL*4 HX0(2,2,2), HX1(2,2), HX2(2), HY0(2,2,2), HY1(2,2), + HY2(2), HZ0(2,2,2), HZ1(2,2), HZ2(2), HX_TMP, FACF INTEGER N_RHO, N_PHI, I, J, K, N_Z DATA PI /3.1415927/ DATA FAC /57.295779/ C Field deflects pos. charged particles to C larger (smaller) polar angles: FACF=+10 ( FACT=-10) DATA FACF / -10.0 / C IF(IMAP.EQ.0) GOTO 1 ! field forced to 0 **** XV=3-dim vector X_B = XV(1)/100.0 Y_B = XV(2)/100.0 Z_B = XV(3)/100.0 C rho = sqrt(x_b*x_b + y_b*y_b) IF(rho.eq.0.0) goto 1 IF((z_b.LE.z_min_field1) .OR. 1 (Z_B.GE.z_max_FIELD1-DELTA_Z_FIELD1/2.)) GOTO 1 IF(rho.GE.rho_max_FIELD1-DELTA_RHO_FIELD1/2.) GOTO 1 n_z = INT((z_b-z_min_field1)/delta_z_field1) n_rho = INT((rho-rho_min_field1)/delta_rho_field1) *** Phi ranging from 0-180 deg (GEANT coordinates) cphi = x_b/rho sphi = y_b/rho phi0 = acos(cphi) IF(y_b.LT.0.0) phi0 = 2*pi - phi0 phi = phi0*fac phi_gap = mod(phi+delta_phi_gap1,delta_phi_gap1)-delta_phi2_gap1 phi_gap_abs = abs(phi_gap) n_phi = phi_gap_abs/delta_phi_field1 weight_z = (z_b-n_z*delta_z_field1-z_min_field1)/delta_z_field1 weight_rho = (rho - n_rho*delta_rho_field1 - rho_min_field1)/ & delta_rho_field1 IF(n_phi.GE.n_phi_max) n_phi=n_phi_max-1 weight_phi = & (phi_gap_abs - n_phi*delta_phi_field1)/delta_phi_field1 C do 4 i=1,2 do 3 j=1,2 do 2 k=1,2 Hx0(i,j,k) = fx1(n_phi+i,n_z+j,n_rho+k) Hy0(i,j,k) = fy1(n_phi+i,n_z+j,n_rho+k) Hz0(i,j,k) = fz1(n_phi+i,n_z+j,n_rho+k) 2 CONTINUE Hx1(i,j) = (Hx0(i,j,2)-Hx0(i,j,1))*weight_rho + Hx0(i,j,1) Hy1(i,j) = (Hy0(i,j,2)-Hy0(i,j,1))*weight_rho + Hy0(i,j,1) Hz1(i,j) = (Hz0(i,j,2)-Hz0(i,j,1))*weight_rho + Hz0(i,j,1) 3 CONTINUE Hx2(i) = (Hx1(i,2)-Hx1(i,1))*weight_z + Hx1(i,1) Hy2(i) = (Hy1(i,2)-Hy1(i,1))*weight_z + Hy1(i,1) Hz2(i) = (Hz1(i,2)-Hz1(i,1))*weight_z + Hz1(i,1) 4 CONTINUE C Hx_tmp = (Hx2(2)-Hx2(1))*weight_phi + Hx2(1) Hy = (Hy2(2)-Hy2(1))*weight_phi + Hy2(1) Hz = (Hz2(2)-Hz2(1))*weight_phi + Hz2(1) C IF(phi_gap.LT.0.0) then Hx_tmp=-Hx_tmp Hz=-Hz ENDIF C /* rotate coordinate system by phi */ Hx = Hx_tmp*cphi - Hy*sphi Hy = Hx_tmp*sphi + Hy*cphi C 6 CONTINUE HV(1) = HX*FACF*FPOL HV(2) = Hy*FACF*FPOL HV(3) = HZ*FACF*FPOL 1 Hx = 0.0 Hy = 0.0 Hz = 0.0 return end