1 ! { dg-options "-O3 -fgraphite-identity -floop-nest-optimize " }
9 public
:: mutual_ind_quad_cir_coil
11 integer, parameter, private
:: longreal
= selected_real_kind(15,90)
12 real (kind
= longreal
), parameter, private
:: pi
= 3.141592653589793_longreal
13 real (kind
= longreal
), parameter, private
:: small
= 1.0e-10_longreal
17 subroutine mutual_ind_quad_cir_coil (r_coil
, x_coil
, y_coil
, z_coil
, h_coil
, n_coil
, &
18 rotate_coil
, m
, mu
, l12
)
19 real (kind
= longreal
), intent(in
) :: r_coil
, x_coil
, y_coil
, z_coil
, h_coil
, n_coil
, &
21 real (kind
= longreal
), dimension(:,:), intent(in
) :: rotate_coil
22 integer, intent(in
) :: m
23 real (kind
= longreal
), intent(out
) :: l12
24 real (kind
= longreal
), dimension(3,3) :: rotate_quad
25 real (kind
= longreal
), dimension(9), save :: x2gauss
, y2gauss
, w2gauss
, z1gauss
, &
27 real (kind
= longreal
) :: xxvec
, xyvec
, xzvec
, yxvec
, yyvec
, yzvec
, zxvec
, zyvec
, &
28 zzvec
, magnitude
, l12_lower
, l12_upper
, dx
, dy
, dz
, theta
, &
29 a
, b1
, b2
, numerator
, denominator
, coefficient
, angle
30 real (kind
= longreal
), dimension(3) :: c_vector
, q_vector
, rot_c_vector
, &
31 rot_q_vector
, current_vector
, &
32 coil_current_vec
, coil_tmp_vector
34 logical, save :: first
= .true
.
37 theta
= pi
*real(i
,longreal
)/real(m
,longreal
)
38 c_vector(1) = r_coil
* cos(theta
)
39 c_vector(2) = r_coil
* sin(theta
)
40 coil_tmp_vector(1) = -sin(theta
)
41 coil_tmp_vector(2) = cos(theta
)
42 coil_tmp_vector(3) = 0.0_longreal
43 coil_current_vec(1) = dot_product(rotate_coil(1,:),coil_tmp_vector(:))
44 coil_current_vec(2) = dot_product(rotate_coil(2,:),coil_tmp_vector(:))
45 coil_current_vec(3) = dot_product(rotate_coil(3,:),coil_tmp_vector(:))
47 c_vector(3) = 0.5 * h_coil
* z1gauss(j
)
48 rot_c_vector(1) = dot_product(rotate_coil(1,:),c_vector(:)) + dx
49 rot_c_vector(2) = dot_product(rotate_coil(2,:),c_vector(:)) + dy
50 rot_c_vector(3) = dot_product(rotate_coil(3,:),c_vector(:)) + dz
52 q_vector(1) = 0.5_longreal
* a
* (x2gauss(k
) + 1.0_longreal
)
53 q_vector(2) = 0.5_longreal
* b1
* (y2gauss(k
) - 1.0_longreal
)
54 q_vector(3) = 0.0_longreal
55 rot_q_vector(1) = dot_product(rotate_quad(1,:),q_vector(:))
56 rot_q_vector(2) = dot_product(rotate_quad(2,:),q_vector(:))
57 rot_q_vector(3) = dot_product(rotate_quad(3,:),q_vector(:))
58 numerator
= w1gauss(j
) * w2gauss(k
) * &
59 dot_product(coil_current_vec
,current_vector
)
60 denominator
= sqrt(dot_product(rot_c_vector
-rot_q_vector
, &
61 rot_c_vector
-rot_q_vector
))
62 l12_lower
= l12_lower
+ numerator
/denominator
66 l12
= coefficient
* (b1
* l12_lower
+ b2
* l12_upper
)
67 end subroutine mutual_ind_quad_cir_coil