build.py

import numpy as np
import logging

logger = logging.getLogger(__name__)
class BuildModule:


    def build(self, model):
        self.calculate_radial_build(model)
        self.calculate_vertical_build(model)
        (
            model.build_data.radius_beam_tangency,
            model.build_data.radius_beam_tangency_max,
        ) = self.calculate_beam_port_size(
            f_r_beam_tangency_rmajor = model.build_data.f_r_beam_tangency_rmajor,
            rmajor = model.geom_data.rmajor,
            n_tf_coil = model.tf_coil_data.n_tf_coil,
            dx_tf_outboard_out = model.tf_coil_data.dx_tf_outboard_out,
            dr_tf_outboard = model.build_data.dr_tf_outboard,
            dx_beam_duct = model.build_data.dx_beam_duct,
            dx_beam_shield = model.build_data.dx_beam_shield,
            r_tf_outboard_mid = model.build_data.r_tf_outboard_mid,
        )
 
    def calculate_radial_build(self, model):
        
        """
        本方法用于确定装置的径向结构构建(radial build)。
        它负责计算装置各部位的厚度、半径、面积等关键径向参数。
        如果 output 标志为 True,可以输出计算结果。

        参数:
            output (bool): 是否输出计算结果

        返回:
            None

        """

        # 在中平面上看
        # inboard part: R < R0
        # outboard part: R > R0
        # 变量命名规则：
        # r_<part name>
        # dr_<part name>_<direction>: 厚度，direction: inboard, outboard
        # dr_gap_<part name1>_<part name2>_<direction>: 两个部分之间的间隙，direction: inboard, outboard
        # R = 0：装置的中心轴，Machine Centreline
        # 装置正中心的空隙，Machine bore
        # 中心螺线管，CS
        # 中心螺线管预紧结构
        # gap
        # TF 线圈内侧腿
        # gap
        # 真空室及屏蔽层(屏蔽层可能在真空室的双层壁之间，或在真空室外部)（热屏蔽层+gap+真空室和（中子）屏蔽层）
            # dr_shld_thermal_inboard: 内侧热屏蔽层厚度
            # dr_gap_shld_thermal_vv_inboard: 热屏蔽层-真空室间隙
            # dr_vv_inboard: 内侧真空室厚度
            # dr_shield_inboard: 内侧屏蔽层厚度
            #
            # shld_thermal: 热屏蔽层, shield: 屏蔽层(和真空室紧贴)
        # gap(作用：热隔离与热膨胀，远程维护与拆装，管路与传感器走线，电磁力缓冲)
        # 内侧包层
        # 内侧第一壁
        # 刮削层（第一壁和等离子体的gap）
        # 等离子体
        # 刮削层
        # 外侧第一壁
        # 外侧包层
        # gap
        # 外侧真空室及屏蔽层
        # gap        
        # TF 线圈外侧支腿
        # gap
        # 杜瓦/低温恒温器


        # ------------------------------------------------ 中心螺线管预紧结构厚度计算 ------------------------------------------------
        # 把预紧结构看成包在 CS 上的薄壳（内外都有预紧壳），只考虑环向（hoop）承载：
        # 环向面积：A = 2 * pi * (r_cs_precomp_inboard + r_cs_precomp_outboard) * dr_cs_precomp
        # f_area_cs_precomp: 预紧结构承载因子，通常取 0.5-0.8, 实际承载面积占CS内外表面积的比例
        # 实际承载面积 = A * f_area_cs_precomp
        # need: i_cs_precomp（是否启用预紧结构）, i_tf_inside_cs（TF线圈是否在CS内部）
        #       force_cs_precomp（预紧力）, sigma_allow_precomp（结构许用应力）, f_area_cs_precomp（预紧结构承载面积因子）
        #       厚度: dr_bore, dr_cs, 
        #       maybe need: dr_tf_inboard, dr_gap_cs_tf_inboard 

        if model.build_data.i_cs_precomp == 1:
            area = model.build_data.force_cs_precomp / (model.build_data.sigma_allow_precomp * model.build_data.f_area_cs_precomp)
            if model.build_data.i_tf_inside_cs == 0:
                # TF 线圈在 cs 外部，常见布置
                # 计算的是最小的厚度
                r_cs_precomp = 2 * model.build_data.dr_bore + model.build_data.dr_cs
            else:
                # TF 线圈在 cs 内部
                r_cs_precomp = (2* (model.build_data.dr_bore + model.build_data.dr_tf_inboard + model.build_data.dr_gap_cs_tf) + 
                    model.build_data.dr_cs)                
            model.build_data.dr_cs_precomp = area / (2.0e0 * np.pi * r_cs_precomp)
        else:
            model.build_data.dr_cs_precomp = 0.0e0

        # ------------------------------------------------ blanket 厚度计算 ------------------------------------------------

        if model.fwbs_data.blkt_model != 0:
            # 分层 blanket 模型
            # blanket = 增值区（breeding unit） + 盒体（manifold） + 后衬板（back plate）
            # 增值区：用于补偿热膨胀，通常比盒体和底板薄
            # 盒体：用于支撑增值区，通常比增值区厚
            # 底板：用于支撑盒体，通常比盒体厚
            # 盒体和底板通常由钨合金制成，具有高熔点、高强度、高导电性等特性
            model.build_data.dr_blkt_inboard = (
                model.build_data.dr_blkt_breeding_inboard + model.build_data.dr_blkt_manifold_inboard + model.build_data.dr_blkt_backplate_inboard)
            model.build_data.dr_blkt_outboard = (
                model.build_data.dr_blkt_breeding_outboard + model.build_data.dr_blkt_manifold_outboard + model.build_data.dr_blkt_backplate_outboard)
            # 考虑一下这里关于 dz 的是不是可以放到垂直方向的计算里？
            model.build_data.dz_shield_upper = 0.5e0 * (
                model.build_data.dr_shield_inboard + model.build_data.dr_shield_outboard)
        
        model.build_data.dz_blkt_upper = 0.5e0 * (
            model.build_data.dr_blkt_inboard + model.build_data.dr_blkt_outboard)
        # model.build_data.dz_blkt_lower = model.build_data.dz_blkt_upper

        """ 不用分层模型，包层厚度如何计算？"""


        """# 若为单 null 配置（如下单分离器），确保前壁与等离子体间隙 dz_fw_plasma_gap 不小于内外侧间隙平均 (?)"""
        if model.prof_data.i_single_null == 1:
            model.build_data.dz_gap_fw_plasma = max(model.build_data.dz_gap_fw_plasma,
                0.5e0 *(model.build_data.dr_gap_fw_plasma_inboard + model.build_data.dr_gap_fw_plasma_outboard))
        


        # ------------------------------------------------ TF coil ------------------------------------------------
        # ------------------------------------------------ TF coil 绝对半径和线包厚度计算 ------------------------------------------------
        # # Issue #514: 内侧 TF 线圈肢厚度的特殊处理
        # # 若 tfcoil_variables.dr_tf_wp_with_insulation 被作为迭代变量（编号140），则用该变量和壳体两部分求 dr_tf_inboard
        # if 140 in numerics.ixc[0 : numerics.nvar]:
        #     build_variables.dr_tf_inboard = (
        #         tfcoil_variables.dr_tf_wp_with_insulation
        #         + tfcoil_variables.dr_tf_plasma_case
        #         + tfcoil_variables.dr_tf_nose_case
        #     )
        
        # TF 内侧支腿的内径/外径/平均半径
        model.build_data.r_tf_inboard_in = model.build_data.dr_bore
        if model.build_data.i_tf_inside_cs == 0:
            model.build_data.r_tf_inboard_in += model.build_data.dr_cs + 2*model.build_data.dr_cs_precomp + model.build_data.dr_gap_cs_tf
        # process 里应该只算了一个 dr_cs_precomp ？

        model.build_data.r_tf_inboard_mid = model.build_data.r_tf_inboard_in + 0.5e0 * model.build_data.dr_tf_inboard
        model.build_data.r_tf_inboard_out = model.build_data.r_tf_inboard_in + model.build_data.dr_tf_inboard

        # 线包（WP）含绝缘层的径向厚度 [m]
        """# 仅在不是作为迭代变量（编号140）的时候计算，变量140可能由外部输入/优化""" # 但这里没有考虑这点
        # if 140 not in numerics.ixc[0 : numerics.nvar]:
        model.build_data.dr_tf_wp_with_insulation = (model.build_data.dr_tf_inboard 
            - model.build_data.dr_tf_inboard_plasma_case # 去掉内侧支腿面向等离子体的壳体厚度
            - model.build_data.dr_tf_inboard_nose_case) # 去掉内侧支腿鼻端壳体厚度（即中心柱一侧的厚度）


        model.tf_coil_data.r_tf_wp_inboard_in = model.build_data.r_tf_inboard_in + model.build_data.dr_tf_inboard_nose_case
        model.tf_coil_data.r_tf_wp_inboard_out = model.build_data.r_tf_inboard_out - model.build_data.dr_tf_inboard_plasma_case


        # ------------------------------------------------ 中心柱（centre post） 顶部半径计算 ------------------------------------------------
        # 中心柱是从中心到 TF 线圈 WP 绕组的一系列结构的总和
        """ 这里是否也要考虑 TF 线圈和 CS 的内外关系？"""
        # 中心柱顶部半径： r_cp_top
        r_cp_top = (model.geom_data.rmajor - model.geom_data.rminor * model.geom_data.delta # 等离子体顶部大半径
            - model.build_data.dr_gap_fw_plasma_inboard # 刮削层（第一壁和等离子体的gap）
            - model.build_data.dr_fw_inboard # 内侧第一壁厚度
            - model.build_data.dr_blkt_inboard # 内侧包层厚度
            - model.build_data.dr_gap_shield_blkt_inboard # 屏蔽层-包层间隙
            - model.build_data.dr_shield_inboard # 内侧屏蔽层厚度
            - model.build_data.dr_vv_inboard # 内侧真空室厚度
            - model.build_data.dr_gap_shld_thermal_vv_inboard # 热屏蔽层-真空室间隙
            - model.build_data.dr_shld_thermal_inboard # 内侧热屏蔽层厚度
            - model.build_data.dr_gap_tf_shld_thermal_inboard # TF 线圈-热屏蔽层间隙
            - model.build_data.dr_tf_inboard_plasma_case)+model.build_data.dr_tf_top
        # PROCESS 没有算 内侧真空室厚度 和 热屏蔽层-真空室间隙，但第一壁和等离子体间隙 * 3
        # model.build_data.dr_tf_top 对应 PROCESS 的 tfcoil_variables.drtop

        if model.geom_data.itart == 1 and model.tf_coil_data.i_tf_superconductor != 1: # 低环径比 + 非超导TF线圈
            if model.build_data.i_r_cp_top == 0:
                # 采用等离子体截面形状的方式设置中心柱顶端半径
                # TF 线圈内侧腿 + gap + 真空室及屏蔽层 + gap + 内侧包层 + gap + 内侧第一壁 + 刮削层 + 等离子体
                model.build_data.r_cp_top = r_cp_top
            elif model.build_data.i_r_cp_top == 1:
                # 作为输入
                pass

            elif model.build_data.i_r_cp_top == 2:
                # 按中面外半径的比例计算
                model.build_data.r_cp_top = (
                    model.build_data.f_r_cp_top * model.build_data.r_tf_inboard_out
                )

            # 如果计算得到的 r_cp_top 太小，会强行修正到 TF线圈外径的1.01倍，并发出警告
            if model.build_data.r_cp_top < 1.01e0 * model.build_data.r_tf_inboard_out:
                logger.error(
                    "TF CP top radius (r_cp_top) replaced by 1.01*r_tf_inboard_out -> potential top rbuild issue"
                    f"{model.build_data.r_cp_top=} {model.build_data.r_tf_inboard_out=}"
                )

                # 强行修正 r_cp_top
                model.build_data.r_cp_top = model.build_data.r_tf_inboard_out * 1.01e0

            if model.build_data.i_r_cp_top != 2:
            # 计算顶部和中面 TF线圈中心柱半径的比值
                model.build_data.f_r_cp_top = model.build_data.r_cp_top / model.build_data.r_tf_inboard_out
        else:  # 以上条件不满足，直接用 r_tf_inboard_out
            model.build_data.r_cp_top = model.build_data.r_tf_inboard_out
            model.build_data.f_r_cp_top = 1.0e0


        # 检查顶部半径是否大于上面用等离子体形状计算得到的最大值，过大会error
        if model.build_data.i_r_cp_top != 0 and (model.build_data.r_cp_top > r_cp_top):
            logger.error(
                "TF CP top radius (r_cp_top) larger that its value determined with plasma shape"
                f"{model.build_data.r_cp_top=} {r_cp_top=}")


        if model.geom_data.itart == 1:
            model.tf_coil_data.dr_tf_wp_top = model.build_data.r_cp_top - model.build_data.r_tf_inboard_in - model.tf_coil_data.dr_tf_nose_case_inboard - model.tf_coil_data.dr_tf_plasma_case_inboard


        # ------------------------------------------------ 内侧真空室计算 ------------------------------------------------
        if model.build_data.i_tf_inside_cs == 0:
            model.build_data.r_vv_inboard_out = (model.build_data.r_tf_inboard_out + model.build_data.dr_gap_tf_shld_thermal_inboard 
                + model.build_data.dr_shld_thermal_inboard + model.build_data.dr_gap_shld_thermal_vv_inboard + model.build_data.dr_vv_inboard)
        else:
            model.build_data.r_vv_inboard_out = (model.build_data.r_tf_inboard_out 
                + model.build_data.dr_gas_cs_tf + model.build_data.dr_cs + 2* model.build_data.dr_cs_precomp + model.build_data.dr_gap_cs_shld_thermal_inboard 
                + model.build_data.dr_shld_thermal_inboard + model.build_data.dr_gap_shld_thermal_vv_inboard + model.build_data.dr_vv_inboard)
                # process 里 dr_cs_precomp 没有 *2，dr_cs_shld_gap 没有，而是有 dr_gap_tf_shld_thermal_inboard

        # 内侧中子屏蔽层的内径/外径
        model.build_data.r_shield_inboard_in = model.build_data.r_vv_inboard_out
        model.build_data.r_shield_inboard_out = model.build_data.r_shield_inboard_in + model.build_data.dr_shield_inboard

        # 等离子体中心半径（rmajor）
        model.build_data.rmajor = (model.build_data.r_shield_inboard_out
            + model.build_data.dr_gap_shield_blkt_inboard + model.build_data.dr_blkt_inboard
            + model.build_data.dr_fw_inboard + model.build_data.dr_gap_fw_plasma_inboard
            + model.geom_data.rminor)
        # 需要约束

        # 内屏蔽层内边缘半径 = 主半径 - 其他
        # 这个变量计算的意义是什么？约束可以由rmajor保证，
        # model.build_data.r_shield_inboard_in = (model.build_data.rmajor - model.geom_data.rminor
        #     - model.build_data.dr_gap_fw_plasma_inboard - model.build_data.dr_fw_inboard - model.build_data.dr_blkt_inboard
        #     - model.build_data.dr_gap_shield_blkt_inboard - model.build_data.dr_shield_inboard)
            # process 里没有算 dr_gap_shield_blkt

        # 外屏蔽层外边缘半径 = 主半径 + 其他
        model.build_data.r_shield_outboard_out = (model.build_data.rmajor + model.geom_data.rminor
            + model.build_data.dr_gap_fw_plasma_outboard + model.build_data.dr_fw_outboard + model.build_data.dr_blkt_outboard 
            + model.build_data.dr_gap_shield_blkt_outboard + model.build_data.dr_shield_outboard)
            # process 里没有算 dr_gap_shield_blkt
       
        
        # ------------------------------------------------ TF 线圈厚度计算 ------------------------------------------------
        # 外侧 TF 线圈肢厚度 = 内侧 TF 线圈肢厚度 * 因子（超导支撑时取1）
        if model.tf_coil_data.i_tf_superconductor != 1:
            model.build_data.dr_tf_outboard = (
                model.build_data.f_dr_tf_outboard_inboard * model.build_data.dr_tf_inboard
            )
        else:
            model.build_data.dr_tf_outboard = model.build_data.dr_tf_inboard


        # 外侧TF线圈肢中心半径（外屏蔽层外径 + 其他）
        model.build_data.r_tf_outboard_mid = (
            model.build_data.r_shield_outboard_out + model.build_data.dr_vv_outboard + model.build_data.dr_gap_shld_thermal_vv_outboard + model.build_data.dr_shld_thermal_outboard
            + model.build_data.dr_gap_tf_shld_thermal_outboard + 0.5e0 * model.build_data.dr_tf_outboard)
        # process 里用 gapomin 替换了 model.build_data.dr_gap_vv_shld_outboard ，多一个 dr_shld_blkt_gap
        model.build_data.r_tf_outboard_in = model.build_data.r_tf_outboard_mid - 0.5e0 * model.build_data.dr_tf_outboard
        model.build_data.r_tf_outboard_out = model.build_data.r_tf_outboard_mid + 0.5e0 * model.build_data.dr_tf_outboard


        # TF线圈可用水平口空间，外TF线圈肢半径-内TF线圈肢半径
        model.build_data.dr_tf_inner_bore = (
            model.build_data.r_tf_outboard_mid - 0.5e0 * model.build_data.dr_tf_outboard
        ) - (model.build_data.r_tf_inboard_mid - 0.5e0 * model.build_data.dr_tf_inboard)
        # 需要核实这里具体是用 外支内径 - 内支内径（process），还是 外支内径 - 内支外径

        # ------------------------------------------------ 等离子体出射面法向的TF线圈磁场 ripple ------------------------------------------------
        # 计算等离子体出射面法向的TF线圈磁场 ripple，对ripple过大的情况自动修正TF线圈外半径
        (delta, r_tf_outboard_mid_min, flag) = self.calculate_toroidal_ripple(
            i_tf_superconductor=model.tf_coil_data.i_tf_superconductor,
            i_tf_wp_geom=model.tf_coil_data.i_tf_wp_geom,
            r_tf_wp_inboard_in=model.tf_coil_data.r_tf_wp_inboard_in,
            r_tf_wp_inboard_mid=model.tf_coil_data.r_tf_wp_inboard_mid,
            r_tf_wp_inboard_out=model.tf_coil_data.r_tf_wp_inboard_out,
            dx_tf_wp_primary_toroidal=model.tf_coil_data.dx_tf_wp_primary_toroidal,
            dx_tf_wp_insulation=model.tf_coil_data.dx_tf_wp_insulation,
            dx_tf_wp_insertion_gap=model.tf_coil_data.dx_tf_wp_insertion_gap,
            r_tf_inboard_out=model.build_data.r_tf_inboard_out,
            r_tf_outboard_mid=model.build_data.r_tf_outboard_mid,
            n_tf_coil=model.tf_coil_data.n_tf_coil,
            i_tf_shape=model.tf_coil_data.i_tf_shape,
            rmajor=model.geom_data.rmajor,
            rminor=model.geom_data.rminor,
            ripple_b_tf_plasma_edge_max=model.tf_coil_data.ripple_b_tf_plasma_edge_max)


        # 若ripple约束要求的外侧半径更大，则自动扩展结构并重新算空间
        if r_tf_outboard_mid_min > model.build_data.r_tf_outboard_mid:
            model.build_data.r_tf_outboard_mid = r_tf_outboard_mid_min
            # 修正真空室与外屏蔽层之间的间隙（重算以配合新的TF半径）
            model.build_data.dr_gap_shld_thermal_vv_outboard = (
                model.build_data.r_tf_outboard_mid
                - 0.5e0 * model.build_data.dr_tf_outboard
                - model.build_data.dr_gap_tf_shld_thermal_outboard
                - model.build_data.dr_shld_thermal_outboard
                - model.build_data.dr_gap_shld_thermal_vv_outboard # process 这里是用的 dr_shld_blkt_gap ,其余一样
                - model.build_data.dr_vv_outboard
                - model.build_data.r_shield_outboard_out)
            # 上面的计算里只有 r_tf_outboard_mid 和 dr_tf_inner_bore 依赖于 dr_gap_shld_thermal_vv_outboard 
        
            model.build_data.dr_tf_inner_bore = (
                model.build_data.r_tf_outboard_mid - 0.5e0 * model.build_data.dr_tf_outboard
            ) - (model.build_data.r_tf_inboard_mid - 0.5e0 * model.build_data.dr_tf_inboard)
            # 和上面计算一样的问题，需确定是外支内径 - 内支内径，还是 外支内径 - 内支外径

        # else:
        #     # 如果不用扩展，恢复到默认最小间隙
        #     build_variables.dr_shld_vv_gap_outboard = build_variables.gapomin
        # 如果不用扩展，dr_gap_shld_thermal_vv_outboard 保持不变，这里删除了 gapomin ，直接用 dr_gap_shld_thermal_vv_outboard 存储最小间隙，

        # 再次修正ripple（假如参数变化）
        (model.tf_coil_data.ripple_b_tf_plasma_edge, r_tf_outboard_mid_min, model.tf_coil_data.ripflag) = self.calculate_toroidal_ripple(
            i_tf_superconductor=model.tf_coil_data.i_tf_superconductor,
            i_tf_wp_geom=model.tf_coil_data.i_tf_wp_geom,
            r_tf_wp_inboard_in=model.tf_coil_data.r_tf_wp_inboard_in,
            r_tf_wp_inboard_mid=model.tf_coil_data.r_tf_wp_inboard_mid,
            r_tf_wp_inboard_out=model.tf_coil_data.r_tf_wp_inboard_out,
            dx_tf_wp_primary_toroidal=model.tf_coil_data.dx_tf_wp_primary_toroidal,
            dx_tf_wp_insulation=model.tf_coil_data.dx_tf_wp_insulation,
            dx_tf_wp_insertion_gap=model.tf_coil_data.dx_tf_wp_insertion_gap,
            r_tf_inboard_out=model.build_data.r_tf_inboard_out,
            r_tf_outboard_mid=model.build_data.r_tf_outboard_mid,
            n_tf_coil=model.tf_coil_data.n_tf_coil,
            i_tf_shape=model.tf_coil_data.i_tf_shape,
            rmajor=model.geom_data.rmajor,
            rminor=model.geom_data.rminor,
            ripple_b_tf_plasma_edge_max=model.tf_coil_data.ripple_b_tf_plasma_edge_max)

        """
        # ------------------------------------------------ 第一壁面积计算 ------------------------------------------------
        # 这部分没有核实修改
        # 这部分不放在 radial build 里
        # 计算第一壁内表面（即靠近等离子体侧）的“半高度”，用于面积计算 
        hbot = (
            build_variables.z_plasma_xpoint_lower             # 等离子体下X点z坐标
            + build_variables.dz_xpoint_divertor              # X点到底部偏移
            + divertor_variables.dz_divertor                  # 整个偏滤器的高度
            - build_variables.dz_blkt_upper                   # 上部包层厚度
            - 0.5e0 * (build_variables.dr_fw_inboard + build_variables.dr_fw_outboard) # 两侧第一壁平均厚度
        )
        if (
            physics_variables.n_divertors == 2
        ):  # 双偏滤器（即单偏滤器标记为0），顶/底半高相等
            htop = hbot
        else:
            # 单偏滤器，顶部为上X点高度+壁-等离子体间隙
            htop = (
                build_variables.z_plasma_xpoint_upper + build_variables.dz_fw_plasma_gap
            )

        # 计算整体第一壁平均半高
        hfw = 0.5e0 * (htop + hbot)

        # 判断用何种形状计算第一壁面积
        if (physics_variables.itart == 1) or (
            fwbs_variables.i_fw_blkt_vv_shape == 1
        ):  # D形截面，面积用D形壳公式
            # 计算内侧部分的中心主半径到几何外边缘
            r1 = (
                physics_variables.rmajor
                - physics_variables.rminor
                - build_variables.dr_fw_plasma_gap_inboard
            )

            # r2是D形截面内半径到外半径的水平距离
            r2 = (
                physics_variables.rmajor
                + physics_variables.rminor
                + build_variables.dr_fw_plasma_gap_outboard
            ) - r1
            # 利用D形壳公式计算内、外第一壁的满覆盖面积和总面积
            (
                build_variables.a_fw_inboard_full_coverage,
                build_variables.a_fw_outboard_full_coverage,
                build_variables.a_fw_total_full_coverage,
            ) = dshellarea(r1, r2, hfw)

        else:  # 假定截面由两个椭圆描述
            # r1: 椭圆截面中心半径
            r1 = (
                physics_variables.rmajor
                - physics_variables.rminor * physics_variables.triang
            )

            # r2: r1到内侧D形外边缘的距离
            r2 = r1 - (
                physics_variables.rmajor
                - physics_variables.rminor
                - build_variables.dr_fw_plasma_gap_inboard
            )

            # r3: r1到外侧D形外边缘的距离
            r3 = (
                physics_variables.rmajor
                + physics_variables.rminor
                + build_variables.dr_fw_plasma_gap_outboard
            ) - r1

            # 利用椭圆壳面积公式计算两侧总面积
            (
                build_variables.a_fw_inboard_full_coverage,
                build_variables.a_fw_outboard_full_coverage,
                build_variables.a_fw_total_full_coverage,
            ) = eshellarea(r1, r2, r3, hfw)

        #  Apply area coverage factor

        if physics_variables.n_divertors == 2:
            # Double null configuration
            build_variables.a_fw_outboard = (
                build_variables.a_fw_outboard_full_coverage
                * (
                    1.0e0
                    - 2.0e0 * fwbs_variables.f_ster_div_single
                    - fwbs_variables.f_a_fw_outboard_hcd
                )
            )
            build_variables.a_fw_inboard = build_variables.a_fw_inboard_full_coverage * (
                1.0e0 - 2.0e0 * fwbs_variables.f_ster_div_single
            )
        else:
            # Single null configuration
            build_variables.a_fw_outboard = (
                build_variables.a_fw_outboard_full_coverage
                * (
                    1.0e0
                    - fwbs_variables.f_ster_div_single
                    - fwbs_variables.f_a_fw_outboard_hcd
                )
            )
            build_variables.a_fw_inboard = build_variables.a_fw_inboard_full_coverage * (
                1.0e0 - fwbs_variables.f_ster_div_single
            )

        build_variables.a_fw_total = (
            build_variables.a_fw_inboard + build_variables.a_fw_outboard
        )

        if build_variables.a_fw_outboard <= 0.0e0:
            raise ProcessValueError(
                "fhole+f_ster_div_single+f_a_fw_outboard_hcd is too high for a credible outboard wall area",
                f_ster_div_single=fwbs_variables.f_ster_div_single,
                f_a_fw_outboard_hcd=fwbs_variables.f_a_fw_outboard_hcd,
            )

        #

        if output:
            #  Print out device build

            po.oheadr(self.outfile, "Radial Build")

            if build_variables.ripflag != 0:
                po.ocmmnt(
                    self.outfile,
                    "(Ripple result may not be accurate, as the fit was outside",
                )
                po.ocmmnt(self.outfile, " its range of applicability.)")
                po.oblnkl(self.outfile)
                logger.warning(
                    "Ripple result may be inaccurate, as the fit has been extrapolated"
                )

                if build_variables.ripflag == 1:
                    diagnostic = (
                        tfcoil_variables.dx_tf_wp_primary_toroidal
                        * tfcoil_variables.n_tf_coils
                        / physics_variables.rmajor
                    )
                    logger.warning(
                        f"(TF coil ripple calculation) Dimensionless coil width X out of fitted range. {diagnostic=}"
                    )
                elif build_variables.ripflag == 2:
                    logger.warning(
                        f"(TF coil ripple calculation) No of TF coils not between 16 and 20 inclusive {tfcoil_variables.n_tf_coils=}"
                    )
                else:
                    diagnostic = (
                        physics_variables.rmajor + physics_variables.rminor
                    ) / build_variables.r_tf_outboard_mid

                    logger.warning(
                        f"(TF coil ripple calculation) (R+a)/rtot={diagnostic} out of fitted range."
                    )

            po.ovarin(
                self.outfile,
                "TF coil radial placement switch",
                "(i_tf_inside_cs)",
                build_variables.i_tf_inside_cs,
            )
            po.ovarrf(
                self.outfile,
                "Inboard build thickness (m)",
                "(dr_inboard_build)",
                physics_variables.rmajor - physics_variables.rminor,
                "OP ",
            )

            if build_variables.i_tf_inside_cs == 1:
                po.ocmmnt(
                    self.outfile,
                    (
                        "\n (The stated machine dr_bore size is just for the hollow space, "
                    ),
                )
                po.ocmmnt(
                    self.outfile,
                    (
                        "the true dr_bore size used for calculations is dr_bore + dr_tf_inboard + dr_cs_tf_gap)\n"
                    ),
                )
            if (
                build_variables.i_tf_inside_cs == 1
                and tfcoil_variables.i_tf_bucking >= 2
            ):
                po.ocmmnt(
                    self.outfile,
                    "(Bore hollow space has been filled with a solid metal cyclinder to act as wedge support)\n",
                )

            # an array that holds the following information
            # description, variable name, thickness, radius
            radial_build_data = []

            radius = 0.0e0
            radial_build_data.append(["Device centreline", None, 0.0, radius])
            if (
                build_variables.i_tf_inside_cs == 1
                and tfcoil_variables.i_tf_bucking >= 2
            ):
                radius = radius + build_variables.dr_bore

                radial_build_data.append([
                    "Machine dr_bore wedge support cylinder",
                    "dr_bore",
                    build_variables.dr_bore,
                    radius,
                ])
            elif (
                build_variables.i_tf_inside_cs == 1 and tfcoil_variables.i_tf_bucking < 2
            ):
                radius = radius + build_variables.dr_bore
                radial_build_data.append([
                    "Machine dr_bore hole",
                    "dr_bore",
                    build_variables.dr_bore,
                    radius,
                ])
            else:
                radius = radius + build_variables.dr_bore
                radial_build_data.append([
                    "Machine dr_bore",
                    "dr_bore",
                    build_variables.dr_bore,
                    radius,
                ])
            if build_variables.i_tf_inside_cs == 1:
                radius += build_variables.dr_tf_inboard
                radial_build_data.append([
                    "TF coil inboard leg (in dr_bore)",
                    "dr_tf_inboard",
                    build_variables.dr_tf_inboard,
                    radius,
                ])

                radius += build_variables.dr_cs_tf_gap
                radial_build_data.append([
                    "CS precompresion to TF coil radial gap",
                    "dr_cs_tf_gap",
                    build_variables.dr_cs_tf_gap,
                    radius,
                ])

            radius = radius + build_variables.dr_cs
            radial_build_data.append([
                "Central solenoid",
                "dr_cs",
                build_variables.dr_cs,
                radius,
            ])

            radius = radius + build_variables.dr_cs_precomp
            radial_build_data.append([
                "CS precompression",
                "dr_cs_precomp",
                build_variables.dr_cs_precomp,
                radius,
            ])
            if build_variables.i_tf_inside_cs == 0:
                radius = radius + build_variables.dr_cs_tf_gap
                radial_build_data.append([
                    "CS precompresion to TF coil radial gap",
                    "dr_cs_tf_gap",
                    build_variables.dr_cs_tf_gap,
                    radius,
                ])

                radius = radius + build_variables.dr_tf_inboard
                radial_build_data.append([
                    "TF coil inboard leg",
                    "dr_tf_inboard",
                    build_variables.dr_tf_inboard,
                    radius,
                ])

            radius = radius + build_variables.dr_tf_shld_gap
            radial_build_data.append([
                "TF coil inboard leg insulation gap",
                "dr_tf_shld_gap",
                build_variables.dr_tf_shld_gap,
                radius,
            ])

            radius = radius + build_variables.dr_shld_thermal_inboard
            radial_build_data.append([
                "Thermal shield, inboard",
                "dr_shld_thermal_inboard",
                build_variables.dr_shld_thermal_inboard,
                radius,
            ])

            radius = radius + build_variables.dr_shld_vv_gap_inboard
            radial_build_data.append([
                "Thermal shield to vessel radial gap",
                "dr_shld_vv_gap_inboard",
                build_variables.dr_shld_vv_gap_inboard,
                radius,
            ])

            radius += build_variables.dr_vv_inboard
            radial_build_data.append([
                "Inboard vacuum vessel",
                "dr_vv_inboard",
                build_variables.dr_vv_inboard,
                radius,
            ])

            radius += build_variables.dr_shld_inboard
            radial_build_data.append([
                "Inner radiation shield",
                "dr_shld_inboard",
                build_variables.dr_shld_inboard,
                radius,
            ])

            radius = radius + build_variables.dr_shld_blkt_gap
            radial_build_data.append([
                "Gap",
                "dr_shld_blkt_gap",
                build_variables.dr_shld_blkt_gap,
                radius,
            ])

            radius = radius + build_variables.dr_blkt_inboard
            radial_build_data.append([
                "Inboard blanket",
                "dr_blkt_inboard",
                build_variables.dr_blkt_inboard,
                radius,
            ])

            radius = radius + build_variables.dr_fw_inboard
            radial_build_data.append([
                "Inboard first wall",
                "dr_fw_inboard",
                build_variables.dr_fw_inboard,
                radius,
            ])

            radius = radius + build_variables.dr_fw_plasma_gap_inboard
            radial_build_data.append([
                "Inboard scrape-off",
                "dr_fw_plasma_gap_inboard",
                build_variables.dr_fw_plasma_gap_inboard,
                radius,
            ])

            radius = radius + physics_variables.rminor
            radial_build_data.append([
                "Plasma geometric centre",
                "rminor",
                physics_variables.rminor,
                radius,
            ])

            radius = radius + physics_variables.rminor
            radial_build_data.append([
                "Plasma outboard edge",
                "rminor",
                physics_variables.rminor,
                radius,
            ])

            radius = radius + build_variables.dr_fw_plasma_gap_outboard
            radial_build_data.append([
                "Outboard scrape-off",
                "dr_fw_plasma_gap_outboard",
                build_variables.dr_fw_plasma_gap_outboard,
                radius,
            ])

            radius = radius + build_variables.dr_fw_outboard
            radial_build_data.append([
                "Outboard first wall",
                "dr_fw_outboard",
                build_variables.dr_fw_outboard,
                radius,
            ])

            radius = radius + build_variables.dr_blkt_outboard
            radial_build_data.append([
                "Outboard blanket",
                "dr_blkt_outboard",
                build_variables.dr_blkt_outboard,
                radius,
            ])

            radius = radius + build_variables.dr_shld_blkt_gap
            radial_build_data.append([
                "Gap",
                "dr_shld_blkt_gap",
                build_variables.dr_shld_blkt_gap,
                radius,
            ])

            radius += build_variables.dr_shld_outboard
            radial_build_data.append([
                "Outer radiation shield",
                "dr_shld_outboard",
                build_variables.dr_shld_outboard,
                radius,
            ])

            radius += build_variables.dr_vv_outboard
            radial_build_data.append([
                "Outboard vacuum vessel",
                "dr_vv_outboard",
                build_variables.dr_vv_outboard,
                radius,
            ])

            radius = radius + build_variables.dr_shld_vv_gap_outboard
            radial_build_data.append([
                "Vessel to TF gap",
                "dr_shld_vv_gap_outboard",
                build_variables.dr_shld_vv_gap_outboard,
                radius,
            ])

            radius = radius + build_variables.dr_shld_thermal_outboard
            radial_build_data.append([
                "Ouboard thermal shield",
                "dr_shld_thermal_outboard",
                build_variables.dr_shld_thermal_outboard,
                radius,
            ])

            radius = radius + build_variables.dr_tf_shld_gap
            radial_build_data.append([
                "Gap",
                "dr_tf_shld_gap",
                build_variables.dr_tf_shld_gap,
                radius,
            ])

            radius = radius + build_variables.dr_tf_outboard
            radial_build_data.append([
                "TF coil outboard leg",
                "dr_tf_outboard",
                build_variables.dr_tf_outboard,
                radius,
            ])

            for description, variable, thickness, radius in radial_build_data:
                po.obuild(
                    self.outfile,
                    description,
                    thickness,
                    radius,
                    f"({variable})" if variable else "",
                )

            # use manual index to ensure count is contiguous in the event
            # of a `None` variable component
            index = 0
            for description, variable, thickness, radius in radial_build_data:
                if variable is None:
                    continue

                index += 1

                po.ovarre(
                    self.mfile,
                    f"{description} radial thickness (m)",
                    f"({variable})",
                    thickness,
                )

                po.ovarst(
                    self.mfile,
                    f"Radial build component {index}",
                    f"(radial_label({index}))",
                    f'"{variable}"',
                )
                po.ovarre(
                    self.mfile,
                    f"Radial build cumulative radius {index}",
                    f"(radial_cum({index}))",
                    radius,
                )

            if (current_drive_variables.i_hcd_primary in [5, 8]) or (
                current_drive_variables.i_hcd_secondary in [5, 8]
            ):
                po.ovarre(
                    self.mfile,
                    "Width of neutral beam duct where it passes between the TF coils (m)",
                    "(dx_beam_duct)",
                    current_drive_variables.dx_beam_duct,
                )
        """

    def calculate_vertical_build(self, model):
        model.build_data.dz_gap_tf_shld_thermal_lower = model.build_data.dr_gap_tf_shld_thermal_inboard
        model.build_data.dz_tf_upper = model.build_data.dr_tf_inboard
        model.build_data.dz_tf_lower = model.build_data.dr_tf_inboard
        model.build_data.dz_gap_shield_blkt_upper = model.build_data.dr_gap_shield_blkt_inboard

        # X point 
        # 假设上下对称
        model.build_data.z_plasma_xpoint_upper = model.geom_data.kappa * model.geom_data.rminor
        model.build_data.z_plasma_xpoint_lower = - model.geom_data.kappa * model.geom_data.rminor

        # divertor geometry: 计算偏滤器高度
        height_div = self.divertor_geom(
            rmajor = model.geom_data.rmajor, 
            rminor = model.geom_data.rminor,
            delta_lower = model.geom_data.delta,
            kappa_lower = model.geom_data.kappa,
            thetai = model.geom_data.thetai,
            thetao = model.geom_data.thetao,
            itart = model.geom_data.itart,
            dist_i = model.divertor_data.dist_i, # 内侧交点到 X 点的距离
            dist_o = model.divertor_data.dist_o, # 外侧交点到 X 点的距离
            lenth_i = model.divertor_data.lenth_i, # 内侧靶板长度
            lenth_o = model.divertor_data.lenth_o, # 外侧靶板长度
            beta_i = model.divertor_data.beta_i, # 内侧靶板与磁面切线的夹角
            beta_o = model.divertor_data.beta_o, # 外侧靶板与磁面切线的夹角
        )

        # 保证等离子体和偏滤器之间有一定距离
        if model.build_data.dz_gap_xpoint_divertor < 1.0e-5:
            model.build_data.dz_gap_xpoint_divertor = height_div


        # 实际上是 TF 等离子体一侧(inside) 的最高点的 z 坐标
        # 有偏滤器时（下侧）
        model.build_data.z_tf_inside_max = (
            model.build_data.z_plasma_xpoint_upper
            + model.build_data.dz_gap_xpoint_divertor
            + model.build_data.dz_divertor
            + model.build_data.dz_shield_lower
            + model.build_data.dz_vv_lower
            + model.build_data.dz_gap_shld_thermal_vv_lower
            + model.build_data.dz_shld_thermal_lower
            + model.build_data.dz_gap_tf_shld_thermal_lower
        )

        # 没有偏滤器时， TF 外侧(outside) 的最高点的 z 坐标
        model.build_data.z_tf_outside_max = (
            model.build_data.z_plasma_xpoint_upper
            + model.build_data.dz_gap_fw_plasma
            + model.build_data.dz_fw_upper
            + model.build_data.dz_blkt_upper
            + model.build_data.dz_gap_shield_blkt_upper
            + model.build_data.dz_shield_upper
            + model.build_data.dz_vv_upper
            + model.build_data.dz_gap_shld_thermal_vv_upper
            + model.build_data.dz_shld_thermal_upper
            + model.build_data.dz_gap_tf_shld_thermal_upper
            + model.build_data.dz_tf_upper
        )

        model.build_data.z_tf_upper_lower_midplane = (
            - model.build_data.z_tf_inside_max
            + model.build_data.z_tf_outside_max - model.build_data.dz_tf_upper
        ) / 2.0e0
        # 上下对称，双零
        if model.prof_data.i_single_null == 0:
            model.build_data.z_tf_outside_max = (model.build_data.z_tf_inside_max
                + model.build_data.dz_tf_upper)
            model.build_data.z_tf_upper_lower_midplane = 0.0e0


    def calculate_beam_port_size(self,
        f_r_beam_tangency_rmajor,
        rmajor,
        n_tf_coil,
        dx_tf_outboard_out,
        dr_tf_outboard,
        dx_beam_duct,
        dx_beam_shield,
        r_tf_outboard_mid):
        r_beam_tangency = f_r_beam_tangency_rmajor * rmajor

        omega = 2* np.pi / n_tf_coil

        a = 0.5e0 * dx_tf_outboard_out
        try:
            assert a < np.inf
        except AssertionError:
            logger.exception("dx_tf_outboard_out / 2 is inf. Kludging to 1e10.")
            a = 1e10

        b = dr_tf_outboard
        try:
            assert b < np.inf
        except AssertionError:
            logger.exception("dr_tf_outboard is inf. Kludging to 1e10.")
            b = 1e10

        c = dx_beam_duct + 2.0e0 * dx_beam_shield
        # 窗口宽度 = 束流宽度 + 2.0e0 * 束流管道壁厚

        d = r_tf_outboard_mid - 0.5e0 * b
        # d = r_tf_outboard_inner  # 为什么不用这一行？
        try:
            assert d < np.inf
        except AssertionError:
            logger.exception("r_tf_outboard_inner is inf. Kludging to 1e10.")
            d = 1e10

        e = np.sqrt(a**2 + (b + d)**2)
        f = np.sqrt(a**2 + d**2)

        phi = omega - np.arctan(a/d) - np.arcsin(a/e)
        g = np.sqrt(e**2 + f**2 - 2*e*f*np.cos(phi))
        if g <= c:
            logger.error(
                f"Max beam tangency radius set =0 temporarily; change dx_beam_duct. {g=} {c=}"
            )
        h = np.sqrt(g**2 - c**2)
        eps = np.arcsin(e * np.sin(phi) / g) - np.arctan(h / c)
        radius_beam_tangency_max = f * np.cos(eps) - 0.5e0 * c
        return r_beam_tangency, radius_beam_tangency_max

    def calculate_toroidal_ripple(self, 
        i_tf_superconductor,
        i_tf_wp_geom,
        r_tf_wp_inboard_in,
        r_tf_wp_inboard_mid,
        r_tf_wp_inboard_out,
        dx_tf_wp_primary_toroidal,
        dx_tf_wp_insulation,
        dx_tf_wp_insertion_gap,
        r_tf_inboard_out,
        r_tf_outboard_mid,
        n_tf_coil,
        i_tf_shape,
        rmajor,
        rminor,
        ripple_b_tf_plasma_edge_max):
        """
        计算等离子体外缘的TF线圈磁场 ripple (delta) ,
        计算临界 ripple 对应的 r_tf_outboard_mid 的最小值
        
        flag: 标记D形线圈拟合公式的准确性
            0: 公式在拟合范围内
            1: x 超出拟合范围
            2: n_tf_coils 超出拟合范围
            3: y 超出拟合范围
        """

        ### 绕组包（winding pack）在环向方向的有效导体厚度 ###
        # 超导TF线圈
        if i_tf_superconductor == 1:
            
            r_wp_min = r_tf_wp_inboard_in
            if i_tf_wp_geom == 0:
                r_wp_max = r_wp_min
            elif i_tf_wp_geom == 1:
                r_wp_max = r_tf_wp_inboard_mid
            elif i_tf_wp_geom == 2:
                r_wp_max = r_tf_wp_inboard_out
            # 上面这个 r_wp_min 和 r_wp_max 计算的意义在？

            # Calculated maximum toroidal WP toroidal thickness [m]
            dx_tf_wp_max = dx_tf_wp_primary_toroidal - 2.0*(dx_tf_wp_insulation + dx_tf_wp_insertion_gap)
        # 常规TF线圈
        else:
            r_wp_max = r_tf_inboard_out
            dx_tf_wp_max = 2.0e0 * r_wp_max * np.tan(np.pi / n_tf_coil)
        
        flag = 0 # 标记D形线圈拟合公式的准确性

        ### calculate ripple ###
        y = (rmajor + rminor) / r_tf_outboard_mid
        if i_tf_shape == 1:
            # D 形线圈
            x = dx_tf_wp_max / rmajor * n_tf_coil
            
            c1 = 0.875e0 - 0.0557e0 * x
            c2 = 1.617e0 + 0.0832e0 * x
            delta = 100.0e0 * c1 * y ** (n_tf_coil - c2)

            yy = ripple_b_tf_plasma_edge_max * 0.001e0 / c1 # y ** (n-c2) > 0
            try:
                assert yy > 1e-6
            except AssertionError:
                logger.exception("(TF coil ripple calculation) y ** (n_tf_coils-c2) <= 1e-6. Kludging to 1e-6.")
                yy = 1e-6
            
            r_tf_outboard_mid_min = (rmajor + rminor) / yy ** (1.0e0 / (n_tf_coil - c2))

            try:
                assert r_tf_outboard_mid_min < np.inf
            except AssertionError:
                logger.exception("r_tf_outboard_mid_min is inf. Kludging to a large value (rmajor + rminor)*3.")
                r_tf_outboard_mid_min = (rmajor + rminor) * 3.0e0
            
            if (x<0.737e0) or (x>2.95e0):
                flag = 1
            elif n_tf_coil not in [16, 18, 20]: # 这个判据和PROCESS代码不完全相同(<16 or >20)
                flag = 2
            elif (y<0.7) or (y>0.8):
                flag = 3
            
        elif i_tf_shape == 2:
            # 矩形线圈
            delta = 100.0e0 * y ** n_tf_coil
            r_tf_outboard_mid_min = (rmajor + rminor) / (ripple_b_tf_plasma_edge_max *0.001e0) ** (1.0e0 / n_tf_coil)

        return delta, r_tf_outboard_mid_min, flag

    def divertor_geom(self, 
        rmajor, 
        rminor,
        delta_lower, # 下方的三角度
        kappa_lower,
        thetai, # model.geom_data.thetai
        thetao,
        itart,
        dist_i, # 内侧交点到 X 点的距离
        dist_o, # 外侧交点到 X 点的距离
        lenth_i, # 内侧靶板长度
        lenth_o, # 外侧靶板长度
        beta_i, # 内侧靶板与磁面切线的夹角
        beta_o): # 外侧靶板与磁面切线的夹角

        if itart == 1:
            # 紧凑型
            return 1.75e0 * rminor
        

        # 经过 X 点的磁面切线 与 水平方向 的夹角
        # 即内侧/外侧的圆弧半径 和 垂直线 的夹角
        theo = np.pi/2 - thetai
        thei = np.pi/2 - thetao

        # 分界面与偏滤器靶板碰撞点（交点）的坐标
        # 内侧/外侧交点：i/o
        # 同一靶板的左右端点：l/r
        # 通过下 X 点计算
        r_x = rmajor - delta_lower * rminor
        z_x = - kappa_lower * rminor

        # 与靶板碰撞点 inter
        r_inter_i = r_x - dist_i * np.cos(thei)
        z_inter_i = z_x - dist_i * np.sin(thei)
        r_inter_o = r_x + dist_o * np.cos(theo)
        z_inter_o = z_x - dist_o * np.sin(theo)

        # 靶板端点
        r_end_i_l = r_inter_i + lenth_i / 2 * np.cos(thei + beta_i)
        z_end_i_l = z_inter_i + lenth_i / 2 * np.sin(thei + beta_i)

        r_end_i_r = r_inter_i - lenth_i / 2 * np.cos(thei + beta_i)
        z_end_i_r = z_inter_i - lenth_i / 2 * np.sin(thei + beta_i)

        r_end_o_l = r_inter_o + lenth_o / 2 * np.cos(theo + beta_o)
        z_end_o_l = z_inter_o - lenth_o / 2 * np.sin(theo + beta_o)

        r_end_o_r = r_inter_o - lenth_o / 2 * np.cos(theo + beta_o)
        z_end_o_r = z_inter_o + lenth_o / 2 * np.sin(theo + beta_o)

        return max(z_end_i_l, z_end_o_r) - min(z_end_i_r, z_end_o_l)