Coverage for python/gsfit/database_readers/st40_astra_mdsplus/setup_coils.py: 0%
139 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 13:12 +0000
1# mypy: ignore-errors
2# TODO: need to fix mypy errors
4import typing
5from typing import TYPE_CHECKING
7import mdsthin
8import numpy as np
9import numpy.typing as npt
10from gsfit_rs import Coils
11from scipy.constants import mu_0
13from .astra_coils_reader import astra_coils_reader
15if TYPE_CHECKING:
16 from . import DatabaseReader
18# Read coil geometry from coil.dat
19astra_coil_description = astra_coils_reader()
22def dividing_parallograms(coil_dictionary) -> tuple[int, int]:
23 """
24 The PF coils are described as parallograms, and discretized in `ASTRA/SRC/Trecur_2.f`, in the `equil_astra_tepm` GitLab repository.
26 ```
27 SW = SQRT( NDIVA*WSIZE/HSIZE )
28 SH = SQRT( NDIVA*HSIZE/WSIZE )
29 SW = SW + 0.5
30 SH = SH + 0.5
31 C
32 NDIVW = IDINT(SW)
33 NDIVH = IDINT(SH)
34 ```
35 """
36 ndi = coil_dictionary["ndi"]
37 wc = coil_dictionary["wc"]
38 hc = coil_dictionary["hc"]
40 sw = np.sqrt(ndi * wc / hc)
41 sh = np.sqrt(ndi * hc / wc)
42 sw = sw + 0.5
43 sh = sh + 0.5
44 ndivw = int(sw)
45 ndivh = int(sh)
46 return ndivw, ndivh
49def calculate_coil_filament_positions(coil_dictionary) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
50 """
51 Calculate the coil filament positions from the coil description dictionary.
53 :param coil_dictionary: Dictionary containing the coil description
54 :return: Tuple of numpy arrays containing the r and z positions of the coil filaments
56 ```
57 C**************************************************************
58 C SUBROUTINE FOR DIVIDING OF PARALLELOGRAM
59 C ( FOR PFC CROSS-SECTION )
60 C**************************************************************
61 C INPUT DATE:
62 C ----------
63 C RC,ZC - CILINDER COORDINATES OF CENTER OF PARALLELOGRAM
64 C WC - PROJECTION OF THE FIRST SIDE OF PARALLELOGRAM
65 C ON AXIS "R" (IN METER)
66 C HC - PROJECTION OF THE SECOND SIDE OF PARALLELOGRAM
67 C ON AXIS "Z" (IN METER)
68 C AWC - ANGLE BETWEEN THE FIRST SIDE OF PARAL. AND
69 C AXIS "R" (IN DEGREES, IT MUST NOT BE EQUAL 90 ,
70 C IT MUST BE : -90 < AWC < 90 )
71 C AHC - ANGLE BETWEEN THE SECOND SIDE OF PARAL. AND
72 C AXIS "R" (IN DEGREES, IT MUST NOT BE EQUAL 0 ,
73 C IT MUST BE : 0 < AHC < 180 )
74 C CURC - CURRENT OF PARALLELOGRAM (IN MA)
75 C NDIV - APPROXIMATE NUMBER OF CELLS OF DIVIDING:
76 C NDIV=0 - A SPECIAL CASE: AUTOMATICALLY NDIVRE=1,
77 C RS(1)=RC, ZS(1)=ZC, PS=CURC
78 C IF NDIV > 0 THEN WE HAVE THE MOST TOTAL ALGORITHM OF
79 C DIVIDING
80 C IF NDIV < 0 THEN NDIVW AND NDIVH ARE CUT OFF BY ABS(NDIV)
81 C
82 C OUTPUT DATE:
83 C ----------
84 C NDIVRE - REAL NUMBER OF CELLS OF DIVIDING ( = NDIVW*NDIVH )
85 C NDIVW - NUMBER OF DIVIDING OF THE FIRST SIDE OF PARAL.
86 C NDIVH - NUMBER OF DIVIDING OF THE SECOND SIDE OF PARAL.
87 C RS(L),ZS(L) - CILINDER COORDINATES OF CENTERS OF CELLS OF DIVIDING
88 C L = 1,2,...,NDIVRE ! (IN METER)
89 C PS - CURRENT OF EVERY CELL OF DIVIDING (IN MA)
90 C VERS - VERTICAL (OR LINEAR, OR RADIUS) SIZE OF CELL CROSS-S.
91 C HORS - HORIZONTAL SIZE OF CELL CROSS-SECTION
92 C
93 C**************************************************************
94 C
95 SUBROUTINE DIVPAR( RC, ZC, WC, HC, AWC, AHC, CURC, NDIV,
96 * NDIVRE, NDIVW, NDIVH, RS, ZS, PS,
97 * VERS, HORS )
98 C
99 C
100 include 'double.inc'
101 C
102 DIMENSION RS(1), ZS(1)
103 C
104 C**************************************************************
105 SIN(X) = DSIN(X)
106 COS(X) = DCOS(X)
107 ATAN(X) = DATAN(X)
108 SQRT(X) = DSQRT(X)
109 C**************************************************************
110 C
111 IF(NDIV.EQ.0) THEN
112 NDIVW = 1
113 NDIVH = 1
114 NDIVRE = 1
115 RS(1) = RC
116 ZS(1) = ZC
117 PS = CURC
118 VERS = HC
119 HORS = WC
120 RETURN
121 END IF
122 C***************************************
123 C
124 IF(AHC.LT.0) AHC = AHC + 180.
125 C
126 IF(NDIV.GE.0) THEN
127 NDIVA = NDIV
128 ELSE
129 NDIVA = -NDIV
130 END IF
131 C***************************************
132 C PARAMETERS OF PARALLELOGRAM
133 C
134 XX = 1.
135 PI = 4.*ATAN(XX)
136 C
137 AWCR = AWC * PI /180.
138 AHCR = AHC * PI /180.
139 C
140 R0 = RC - 0.5*( WC + HC * COS(AHCR)/SIN(AHCR) )
141 Z0 = ZC - 0.5*( HC + WC * SIN(AWCR)/COS(AWCR) )
142 C
143 WSIZE = WC / COS(AWCR)
144 HSIZE = HC / SIN(AHCR)
145 C
146 WR = WC
147 WZ = WC * SIN(AWCR) / COS(AWCR)
148 HR = HC * COS(AHCR) / SIN(AHCR)
149 HZ = HC
150 C***************************************
151 C CALCULATION NDIVW, NDIVH, NDIVRE, PS
152 C
153 SW = SQRT( NDIVA*WSIZE/HSIZE )
154 SH = SQRT( NDIVA*HSIZE/WSIZE )
155 SW = SW + 0.5
156 SH = SH + 0.5
157 C
158 NDIVW = IDINT(SW)
159 NDIVH = IDINT(SH)
160 C
161 IF(NDIVW.EQ.0) NDIVW = 1
162 IF(NDIVH.EQ.0) NDIVH = 1
163 C
164 IF((NDIV.LT.0).AND.(NDIVW.GT.NDIVA)) NDIVW = NDIVA
165 IF((NDIV.LT.0).AND.(NDIVH.GT.NDIVA)) NDIVH = NDIVA
166 C
167 NDIVRE = NDIVW * NDIVH
168 PS = CURC / NDIVRE
169 C***************************************
170 C CALCULATION RS(L), ZS(L) : L = 1,2,...,NDIVRE
171 C
172 WR = WR / NDIVW
173 WZ = WZ / NDIVW
174 HR = HR / NDIVH
175 HZ = HZ / NDIVH
176 C
177 HORS = WR
178 VERS = HZ
179 C
180 RS(1) = R0 + 0.5*(WR + HR)
181 ZS(1) = Z0 + 0.5*(WZ + HZ)
182 C
183 DO 1 I=1,NDIVW
184 DO 1 J=1,NDIVH
185 L = (I-1)*NDIVH + J
186 RS(L) = RS(1) + (I-1)*WR + (J-1)*HR
187 ZS(L) = ZS(1) + (I-1)*WZ + (J-1)*HZ
188 1 CONTINUE
189 C***************************************
190 C
191 RETURN
192 END
193 ```
194 """
196 # Extract required parameters
197 rc = float(coil_dictionary["rc"]) # center R
198 zc = float(coil_dictionary["zc"]) # center Z
199 wc = float(coil_dictionary["wc"]) # projection of first side on R
200 hc = float(coil_dictionary["hc"]) # projection of second side on Z
201 # Optional angles (degrees). Default to AWC=0, AHC=90 if not provided,
202 # matching the most common cases in the coil table.
203 awc = float(coil_dictionary.get("awc", 0.0))
204 ahc = float(coil_dictionary.get("ahc", 90.0))
205 # Approximate number of cells (NDIV/NDIVA)
206 ndia = int(coil_dictionary.get("ndi", 1))
208 # Handle the special case NDIV == 0
209 if ndia == 0:
210 rs = np.array([rc], dtype=np.float64)
211 zs = np.array([zc], dtype=np.float64)
212 return rs, zs
214 # Convert angles to radians, applying the Fortran adjustment for AHC < 0
215 if ahc < 0.0:
216 ahc = ahc + 180.0
218 pi = np.pi
219 awcr = awc * pi / 180.0
220 ahcr = ahc * pi / 180.0
222 # Guard against invalid angles (avoid division by zero)
223 # In Fortran code: it assumes -90 < AWC < 90 and 0 < AHC < 180
224 # If cos(awcr) or sin(ahcr) are too close to zero, nudge slightly.
225 eps = 1e-12
226 c_aw = np.cos(awcr)
227 s_aw = np.sin(awcr)
228 s_ah = np.sin(ahcr)
229 c_ah = np.cos(ahcr)
230 if abs(c_aw) < eps:
231 c_aw = eps if c_aw >= 0 else -eps
232 if abs(s_ah) < eps:
233 s_ah = eps if s_ah >= 0 else -eps
235 # Fortran mapping
236 r0 = rc - 0.5 * (wc + hc * c_ah / s_ah)
237 z0 = zc - 0.5 * (hc + wc * s_aw / c_aw)
239 wsize = wc / c_aw
240 hsize = hc / s_ah
242 wr = wc
243 wz = wc * s_aw / c_aw
244 hr = hc * c_ah / s_ah
245 hz = hc
247 # Calculate NDIVW, NDIVH per Fortran's DIVPAR (NDIV >= 0 path)
248 sw = np.sqrt(ndia * wsize / hsize) + 0.5
249 sh = np.sqrt(ndia * hsize / wsize) + 0.5
250 ndivw = int(sw)
251 ndivh = int(sh)
252 if ndivw == 0:
253 ndivw = 1
254 if ndivh == 0:
255 ndivh = 1
257 # Real number of divided cells
258 ndivre = ndivw * ndivh
260 # Cell sizes
261 wr_cell = wr / ndivw
262 wz_cell = wz / ndivw
263 hr_cell = hr / ndivh
264 hz_cell = hz / ndivh
266 # First cell center
267 r1 = r0 + 0.5 * (wr_cell + hr_cell)
268 z1 = z0 + 0.5 * (wz_cell + hz_cell)
270 # Generate centers for all cells (Fortran order L = (I-1)*NDIVH + J)
271 rs = np.empty(ndivre, dtype=np.float64)
272 zs = np.empty(ndivre, dtype=np.float64)
273 idx = 0
274 for i in range(1, ndivw + 1):
275 for j in range(1, ndivh + 1):
276 r = r1 + (i - 1) * wr_cell + (j - 1) * hr_cell
277 z = z1 + (i - 1) * wz_cell + (j - 1) * hz_cell
278 rs[idx] = r
279 zs[idx] = z
280 idx += 1
282 return rs, zs
285def setup_coils(
286 self: "DatabaseReader",
287 settings: dict[str, typing.Any],
288 pulseNo: int,
289) -> Coils:
290 """
291 This method initialises the Rust `Coils` class.
293 :param pulseNo: Pulse number, used to read from the database
294 :param settings: Dictionary containing the JSON settings read from the `settings` directory
296 **This method is specific to ST40's ASTRA stored on MDSplus.**
298 See `python/gsfit/database_readers/interface.py` for more details on how a new database_reader should be implemented.
299 """
301 # Initialise the Coils Rust class
302 coils = Coils()
304 # Extract the astra_run_name from settings
305 astra_run_name = settings["GSFIT_code_settings.json"]["database_reader"]["st40_astra_mdsplus"]["workflow"]["astra"]["run_name"]
307 # Connect to MDSplus
308 conn = mdsthin.Connection("smaug")
309 conn.openTree("ASTRA", pulseNo)
311 # ASTRA time
312 time = conn.get(f"\\ASTRA::TOP.{astra_run_name}:TIME").data().astype(np.float64)
314 # PF coil currents
315 # ["IPL", "MC", "PSH", "DIV", "BVL", "BVUT", "BVUB", "CS", "MCVC"]
316 # [ 0, 1, 2, 3, 4, 5, 6, 7 8]
317 currents = conn.get(f"\\ASTRA::TOP.{astra_run_name}.COILS.PSU2PF:I").data().astype(np.float64) * 1.0e6
319 # Add BVLB PF coil
320 coil_name = "BVLB"
321 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
322 # Read in BVL PSU current, and normalise by the number of filaments
323 currents_local = currents[:, 4] * 4.0 * 4.0 / len(coil_r)
324 coils.add_pf_coil(
325 coil_name,
326 coil_r,
327 coil_z,
328 d_r=0.0 * coil_r,
329 d_z=0.0 * coil_z,
330 time=time,
331 measured=currents_local,
332 )
333 # Add BVLT PF coil
334 coil_name = "BVLT"
335 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
336 coils.add_pf_coil(
337 coil_name,
338 coil_r,
339 coil_z,
340 d_r=0.0 * coil_r,
341 d_z=0.0 * coil_z,
342 time=time,
343 measured=currents_local,
344 )
346 # Read in BVUB PSU current, and normalise by the number of filaments
347 coil_name = "BVUB"
348 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
349 currents_local = currents[:, 6] * 4.0 * 6.0 / len(coil_r)
350 coils.add_pf_coil(
351 coil_name,
352 coil_r,
353 coil_z,
354 d_r=0.0 * coil_r,
355 d_z=0.0 * coil_z,
356 time=time,
357 measured=currents_local,
358 )
360 # Read in BVUT PSU current, and normalise by the number of filaments
361 coil_name = "BVUT"
362 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
363 currents_local = currents[:, 5] * 4.0 * 6.0 / len(coil_r)
364 coils.add_pf_coil(
365 coil_name,
366 coil_r,
367 coil_z,
368 d_r=0.0 * coil_r,
369 d_z=0.0 * coil_z,
370 time=time,
371 measured=currents_local,
372 )
374 # Read in SOL PSU current, and normalise by the number of filaments
375 coil_name = "SOL"
376 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
377 currents_local = currents[:, 7] * 2.0 * 95.0 / len(coil_r)
378 coils.add_pf_coil(
379 coil_name,
380 coil_r,
381 coil_z,
382 d_r=0.0 * coil_r,
383 d_z=0.0 * coil_z,
384 time=time,
385 measured=currents_local,
386 )
388 # Add DIVT PF coil, and normalise by the number of filaments
389 coil_name = "DIVT"
390 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
391 currents_local = currents[:, 3] * 4.0 * 7.0 / len(coil_r)
392 coils.add_pf_coil(
393 coil_name,
394 coil_r,
395 coil_z,
396 d_r=0.0 * coil_r,
397 d_z=0.0 * coil_z,
398 time=time,
399 measured=currents_local,
400 )
401 # Add DIVB PF coil
402 coil_name = "DIVB"
403 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
404 coils.add_pf_coil(
405 coil_name,
406 coil_r,
407 coil_z,
408 d_r=0.0 * coil_r,
409 d_z=0.0 * coil_z,
410 time=time,
411 measured=currents_local,
412 )
414 # Add MCT PF coil, and normalise by the number of filaments
415 coil_name = "MCT"
416 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
417 currents_local = -currents[:, 8] * 11.0 / len(coil_r) / 2.0 + currents[:, 1] * 11.0 / len(coil_r)
418 coils.add_pf_coil(
419 coil_name,
420 coil_r,
421 coil_z,
422 d_r=0.0 * coil_r,
423 d_z=0.0 * coil_z,
424 time=time,
425 measured=currents_local,
426 )
427 # Add MCB PF coil
428 coil_name = "MCB"
429 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
430 currents_local = currents[:, 8] * 11.0 / len(coil_r) / 2.0 + currents[:, 1] * 11.0 / len(coil_r)
431 coils.add_pf_coil(
432 coil_name,
433 coil_r,
434 coil_z,
435 d_r=0.0 * coil_r,
436 d_z=0.0 * coil_z,
437 time=time,
438 measured=currents_local,
439 )
441 # Add PSHT PF coil
442 coil_name = "PSHT1"
443 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
444 currents_local = currents[:, 2] * 4.0 / len(coil_r)
445 coils.add_pf_coil(
446 coil_name,
447 coil_r,
448 coil_z,
449 d_r=0.0 * coil_r,
450 d_z=0.0 * coil_z,
451 time=time,
452 measured=currents_local,
453 )
454 # Add PSHB PF coil
455 coil_name = "PSHB1"
456 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
457 coils.add_pf_coil(
458 coil_name,
459 coil_r,
460 coil_z,
461 d_r=0.0 * coil_r,
462 d_z=0.0 * coil_z,
463 time=time,
464 measured=currents_local,
465 )
466 # Add PSHT2 PF coil
467 coil_name = "PSHT2"
468 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
469 coils.add_pf_coil(
470 coil_name,
471 coil_r,
472 coil_z,
473 d_r=0.0 * coil_r,
474 d_z=0.0 * coil_z,
475 time=time,
476 measured=currents_local,
477 )
478 # Add PSHB2 PF coil
479 coil_name = "PSHB2"
480 coil_r, coil_z = calculate_coil_filament_positions(astra_coil_description[coil_name])
481 coils.add_pf_coil(
482 coil_name,
483 coil_r,
484 coil_z,
485 d_r=0.0 * coil_r,
486 d_z=0.0 * coil_z,
487 time=time,
488 measured=currents_local,
489 )
491 # Add TF coil
492 bt_vac = conn.get(f"\\ASTRA::TOP.{astra_run_name}.GLOBAL:BTVAC").data().astype(np.float64) # time-dependent
493 r_reference = 0.5
494 i_rod = bt_vac * (2.0 * np.pi * r_reference) / mu_0 # time-dependent
495 coils.add_tf_coil(
496 time=time,
497 measured=i_rod,
498 )
500 return coils