<!DOCTYPE html>
<html lang="zh-CN">
<head>
    <meta charset="UTF-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <title>Coulomb Barrier Height Calculator (with B<sub>cap</sub>)</title>
    <!-- Chart.js -->
    <script src="https://cdn.jsdelivr.net/npm/chart.js"></script>
    <style>
        body {
            font-family: Arial, sans-serif;
            max-width: 1300px;
            margin: 0 auto;
            padding: 20px;
            background-color: #f5f5f5;
        }
        .container {
            background-color: white;
            border-radius: 8px;
            padding: 20px;
            box-shadow: 0 2px 10px rgba(0,0,0,0.1);
        }
        h1 {
            color: #2c3e50;
            text-align: center;
            margin-bottom: 30px;
        }
        .input-group {
            display: flex;
            flex-wrap: wrap;
            gap: 15px;
            margin-bottom: 20px;
        }
        .input-field {
            flex: 1;
            min-width: 120px;
        }
        label {
            display: block;
            margin-bottom: 5px;
            font-weight: bold;
            color: #34495e;
        }
        input {
            width: 100%;
            padding: 8px;
            border: 1px solid #ddd;
            border-radius: 4px;
            box-sizing: border-box;
        }
        .button-group {
            display: flex;
            gap: 15px;
            margin-top: 10px;
        }
        button, .back-link {
            color: white;
            border: none;
            padding: 10px 20px;
            border-radius: 4px;
            cursor: pointer;
            font-size: 16px;
            flex: 1;
            text-align: center;
            text-decoration: none;
            display: inline-block;
        }
        #calculateBtn {
            background-color: #3498db;
        }
        #calculateBtn:hover {
            background-color: #2980b9;
        }
        .back-link {
            background-color: #4CAF50;
        }
        .back-link:hover {
            background-color: #7f8c8d;
        }
        .results {
            margin-top: 30px;
            display: none;
        }
        .reaction-system {
            text-align: center;
            font-size: 20px;
            font-weight: bold;
            color: #2c3e50;
            margin-bottom: 15px;
            padding: 10px;
            background-color: #f8f9fa;
            border-radius: 5px;
            border-left: 4px solid #3498db;
        }
        .result-table {
            width: 100%;
            border-collapse: collapse;
            margin-top: 15px;
            margin-bottom: 20px;
        }
        .result-table th, .result-table td {
            border: 1px solid #ddd;
            padding: 10px;
            text-align: left;
        }
        .result-table th {
            background-color: #f2f2f2;
        }
        .bcap-table {
            margin-top: 10px;
            margin-bottom: 20px;
        }
        .chart-container {
            margin-top: 30px;
            margin-bottom: 30px;
            position: relative;
            height: 600px;
            width: 100%;
        }
        .note {
            margin-top: 20px;
            font-size: 0.9em;
            color: #7f8c8d;
        }
        .note ul {
            padding-left: 20px;
        }
        .note li {
            margin-bottom: 5px;
        }
        .section-title {
            font-size: 1.2em;
            font-weight: bold;
            margin: 20px 0 10px;
            color: #2c3e50;
            border-bottom: 1px solid #ccc;
            padding-bottom: 5px;
        }
    </style>
</head>
<body>
<div class="container">
    <h1>Coulomb Barrier in Heavy-ion Fusion Reaction v2.1</h1>
    <p>Please input the mass and charge numbers of projectile and target nuclei</p>

    <div class="input-group">
        <div class="input-field">
            <label>A<sub>1</sub></label>
            <input type="number" id="a1" placeholder="e.g., 54" step="1" value="54">
        </div>
        <div class="input-field">
            <label>Z<sub>1</sub></label>
            <input type="number" id="z1" placeholder="e.g., 24" step="1" value="24">
        </div>
        <div class="input-field">
            <label>A<sub>2</sub> </label>
            <input type="number" id="a2" placeholder="e.g., 243" step="1" value="243">
        </div>
        <div class="input-field">
            <label>Z<sub>2</sub></label>
            <input type="number" id="z2" placeholder="e.g., 95" step="1" value="95">
        </div>
    </div>

    <div class="button-group">
        <button id="calculateBtn">Calculate & Plot</button>
        <a href="http://www.imqmd.com/fusion/" class="back-link"> &lt;&lt; <strong> Index </strong> &gt;&gt; </a>
    </div>

    <div class="results" id="results">
        <!-- 反应体系 -->
        <div class="reaction-system" id="reactionSystem">
            Reaction System: <span id="reactionText">-</span>
        </div>

        <!-- 势垒高度表格 (原有) -->
        <table class="result-table">
            <thead>
            <tr>
                <th>Model</th>
                <th>Barrier height (MeV)</th>
                <th>Barrier radius (fm)</th>
                <th>Pocket depth (MeV)</th>
            </tr>
            </thead>
            <tbody>
            <tr>
                <td>Bass80</td>
                <td id="bass80-vb">-</td>
                <td id="bass80-rb">-</td>
                <td id="bass80-bcap">-</td>
           </tr>
            <tr>
                <td>MBass</td>
                <td id="mbass80-vb">-</td>
                <td id="mbass80-rb">-</td>
                <td id="mbass80-bcap">-</td>
            </tr>
            <tr>
                <td>BW91*</td>
                <td id="bw91-vb">-</td>
                <td id="bw91-rb">-</td>
                <td id="bw91-bcap">-</td>
            </tr>
            <tr>
                <td>MWS2</td>
                <td id="mws-vb">-</td>
                <td id="mws-rb">-</td>
                <td id="mws-bcap">-</td>
            </tr>
            </tbody>
        </table>
 

        <!-- 图表容器 -->
        <div class="chart-container">
            <canvas id="potentialChart"></canvas>
        </div>

        <div class="note">
            <p><strong>References: </strong></p>
            <ul>
                <li>Bass80: Bass potential (1980); Lecture Notes in Physics 117 (Berlin: Springer) pp 281</li>
                <li>MBass : Bass80 potential + a repulsive core (Ning Wang 2025-12-9 )</li>
                <li>BW91*: Broglia and Winther potential (1991) + B<sub>cap</sub> correction (Ning Wang 2026-4-15); Nucl. Phys. A <strong>899</strong> (2013) 29;   </li>
                <li>MWS2: Modified Woods-Saxon potential + B<sub>cap</sub> correction; J. L. Tian, N. Wang, Z. X. Li, Chin. Phys. Lett. <strong>24</strong> (2007) 905 [<a href="http://www.imqmd.com/wangning/publication/CPL24-905.pdf">PDF</a>] </li>
            </ul>
        </div>
    </div>
</div>

<script>
    (function() {
        // 元素符号
        const ELEMENT_SYMBOLS = ['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne',
            'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', 'Sc', 'Ti', 'V',
            'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr',
            'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In',
            'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm',
            'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re',
            'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra',
            'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md',
            'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', 'Rg', 'Cn', 'Nh', 'Fl',
            'Mc', 'Lv', 'Ts', 'Og', '119', '120', '121', '122', '123', '124', '125', '126'];

        // 势函数定义 (保持与原代码一致)
        function BW91(A1, Z1, A2, Z2, R) {
            const pi = 3.1415926;
            const e2 = 1.44;
            const gamma0 = 0.95;
            const xks = 1.8;
            const a = 0.63;
            const R1 = 1.233 * Math.pow(A1, 1/3) - 0.98 / Math.pow(A1, 1/3);
            const R2 = 1.233 * Math.pow(A2, 1/3) - 0.98 / Math.pow(A2, 1/3);
            const R12 = R1 * R2 / (R1 + R2);
            const d0 = R1 + R2 + 0.29;
            const gamma = gamma0 * (1 - xks * (A1 - 2 * Z1) / A1 * (A2 - 2 * Z2) / A2);
            const V0 = 16 * pi * gamma * R12 * a;
            const V1 = 1.34 * Math.pow(d0/R, 4) * R12 * R12;
            const VN = (-V0+V1) / (1 + Math.exp((R - d0) / a));
     
			let Vc;
            if (R <= R1 + R2) {
                Vc = 0.5 * Z1 * Z2 * e2 * (3 - Math.pow(R, 2) / Math.pow(R1 + R2, 2)) / (R1 + R2);
            } else {
                Vc = Z1 * Z2 * e2 / R;
            }

			return VN + Vc;
        }

        function Woods_Saxon(A1, Z1, A2, Z2, R) {
            const pi = 3.1415926;
            const gamma0 = 0.95;
            const xks = 1.8;
            const a = 1 / (1.17 * (1 + 0.53 * (1 / Math.pow(A1, 1/3) + 1 / Math.pow(A2, 1/3))));
            const R1 = 1.20 * Math.pow(A1, 1/3) - 0.09;
            const R2 = 1.20 * Math.pow(A2, 1/3) - 0.09;
            const R12 = R1 * R2 / (R1 + R2);
            const d0 = R1 + R2;
            const gamma = gamma0 * (1 - xks * (A1 - 2 * Z1) / A1 * (A2 - 2 * Z2) / A2);
            const V0 = 16 * pi * gamma * R12 * a;
            const VN = -V0 / (1 + Math.exp((R - d0) / a));
            return VN + Z1 * Z2 * 1.44 / R;
        }

        function VBASS80(A1, Z1, A2, Z2, R) {
            const e2 = 1.44;
            const Rs1 = 1.28 * Math.pow(A1, 1/3) - 0.76 + 0.8 / Math.pow(A1, 1/3);
            const Rs2 = 1.28 * Math.pow(A2, 1/3) - 0.76 + 0.8 / Math.pow(A2, 1/3);
            const R1 = Rs1 * (1 - 0.98 / Math.pow(Rs1, 2));
            const R2 = Rs2 * (1 - 0.98 / Math.pow(Rs2, 2));
            const s = R - R1 - R2;
            const R12 = R1 * R2 / (R1 + R2);
            const Vn = R12 / (0.033 * Math.exp(s / 3.5) + 0.007 * Math.exp(s / 0.65));
            const Vc = Z1 * Z2 * e2 / R;
            return -Vn + Vc;
        }

        function MBASS80(A1, Z1, A2, Z2, R) {
            const e2 = 1.44;
            const Rs1 = 1.28 * Math.pow(A1, 1/3) - 0.76 + 0.8 / Math.pow(A1, 1/3);
            const Rs2 = 1.28 * Math.pow(A2, 1/3) - 0.76 + 0.8 / Math.pow(A2, 1/3);
            const R1 = Rs1 * (1 - 0.98 / Math.pow(Rs1, 2));
            const R2 = Rs2 * (1 - 0.98 / Math.pow(Rs2, 2));
            const s = R - R1 - R2;
            const R12 = R1 * R2 / (R1 + R2);
            const Vn = R12 / (0.033 * Math.exp(s / 3.5) + 0.007 * Math.exp(s / 0.65) +
                (0.02 - 0.011 / R) * Math.exp(-R / (R1 + R2)));
            let Vc;
            if (R <= R1 + R2) {
                Vc = 0.5 * Z1 * Z2 * e2 * (3 - Math.pow(R, 2) / Math.pow(R1 + R2, 2)) / (R1 + R2);
            } else {
                Vc = Z1 * Z2 * e2 / R;
            }
            return -Vn + Vc;
        }

        function VMWS(A1, Z1, A2, Z2, R) {
            const r0 = 1.27;
            const c = -1.37;
            const as_param = -44.16;
            const a = 0.75;
            const cs = -0.40;
            const A1_13 = Math.pow(A1, 1/3);
            const A2_13 = Math.pow(A2, 1/3);
            const R12 = (A1_13 * A2_13) / (A1_13 + A2_13);
            const d0 = r0 * (A1_13 + A2_13) + c;
            const V0 = as_param * (1 + cs * ((A1 - 2*Z1)/A1 + (A2 - 2*Z2)/A2)) * R12;
            const V1 = 1.75 * Math.abs(A1-A2) / (A1+A2) * (d0 + c) * (d0 + c) * R12 / (R * R);
            const VN = (V0 + V1) / (1 + Math.exp((R - d0) / a));
            const VC = Z1 * Z2 * 1.44 / R;
            return VN + VC;
        }

        // 计算势垒高度 (保留原有逻辑)
        function calculateBarrierHeights(A1, Z1, A2, Z2) {
            let Vmax_bw91 = -Infinity, Rbw = 0;
            let Vmax_ws = -Infinity, Rws = 0;
            let Vmax_bass80 = -Infinity, Rbass = 0;
            let Vmax_mws = -Infinity, Rfus2 = 0;
            let Vmax_mbass80 = -Infinity, Rmbass = 0;

            const R_start = 1.18 * Math.pow(A1, 1/3) + 1.18 * Math.pow(A2, 1/3);
            const R_end = 20;
            for (let R = R_start; R <= R_end; R += 0.1) {
                const vBw = BW91(A1, Z1, A2, Z2, R);
                if (vBw > Vmax_bw91) { Vmax_bw91 = vBw; Rbw = R; }

                const vWs = Woods_Saxon(A1, Z1, A2, Z2, R);
                if (vWs > Vmax_ws) { Vmax_ws = vWs; Rws = R; }

                const vBass = VBASS80(A1, Z1, A2, Z2, R);
                if (vBass > Vmax_bass80) { Vmax_bass80 = vBass; Rbass = R; }

                const vMws = VMWS(A1, Z1, A2, Z2, R);
                if (vMws > Vmax_mws) { Vmax_mws = vMws; Rfus2 = R; }

                const vMbass = MBASS80(A1, Z1, A2, Z2, R);
                if (vMbass > Vmax_mbass80) { Vmax_mbass80 = vMbass; Rmbass = R; }
            }

            return {
                bw91: { vb: Vmax_bw91, rb: Rbw },
                aw95: { vb: Vmax_ws, rb: Rws },
                bass80: { vb: Vmax_bass80, rb: Rbass },
                mws: { vb: 0.956 * Vmax_mws, rb: Rfus2 },   // MWS2 显示值
                mbass80: { vb: Vmax_mbass80, rb: Rmbass },
                // 同时返回未缩放的MWS最大值位置等，便于后续使用
                rawMwsMax: Vmax_mws,
                rawMwsRb: Rfus2
            };
        }

        // 计算俘获势阱深度 Bcap (以 MWS 的 d0 为界)
        function computeBcap(potentialFunc, A1, Z1, A2, Z2, d0_mws) {
            const step = 0.05;  // 稍精细
            const R_min_global = Math.max(2.5, 0.4 * d0_mws);
            const R_max_global = 22.0;

            // 寻找 R < d0 区域的最小值
            let Vmin_left = Infinity;
            for (let R = R_min_global; R < d0_mws; R += step) {
                const V = potentialFunc(A1, Z1, A2, Z2, R);
                if (V < Vmin_left) Vmin_left = V;
            }

            // 寻找 R > d0 区域的最大值 (势垒)
            let Vmax_right = -Infinity;
            for (let R = d0_mws; R <= R_max_global; R += step) {
                const V = potentialFunc(A1, Z1, A2, Z2, R);
                if (V > Vmax_right) Vmax_right = V;
            }

            // 如果最小值没找到（极少情况），回退处理
            if (Vmin_left === Infinity) Vmin_left = potentialFunc(A1, Z1, A2, Z2, d0_mws - 0.1);
            if (Vmax_right === -Infinity) Vmax_right = potentialFunc(A1, Z1, A2, Z2, d0_mws + 0.1);

            return Vmax_right - Vmin_left;
        }

        // 获取元素符号
        function getElementSymbol(Z) {
            if (Z >= 1 && Z <= ELEMENT_SYMBOLS.length) {
                return ELEMENT_SYMBOLS[Z - 1];
            }
            return `Z=${Z}`;
        }

        // 绘图
        let potentialChart = null;
        function plotPotentialCurves(A1, Z1, A2, Z2, results, reactionText) {
            const R_start = 1.20 * (Math.pow(A1, 1/3) + Math.pow(A2, 1/3)) - 2.9;
            const R_end = 1.20 * (Math.pow(A1, 1/3) + Math.pow(A2, 1/3)) + 5.5;
            const step = 0.1;
            const R_values = [];
            const bw91_values = [];
            const bass80_values = [];
            const mbass80_values = [];

            for (let R = R_start; R <= R_end; R += step) {
                R_values.push(R);
                bw91_values.push(BW91(A1, Z1, A2, Z2, R));
                bass80_values.push(VBASS80(A1, Z1, A2, Z2, R));
                mbass80_values.push(MBASS80(A1, Z1, A2, Z2, R));
            }

            const ctx = document.getElementById('potentialChart').getContext('2d');
            if (potentialChart) potentialChart.destroy();

            potentialChart = new Chart(ctx, {
                type: 'line',
                data: {
                    labels: R_values.map(r => r.toFixed(1)),
                    datasets: [
                        { label: 'BW91*', data: bw91_values, borderColor: '#e74c3c', backgroundColor: 'rgba(231,76,60,0.1)', borderWidth: 2, tension: 0.1, pointRadius: 0 },
                        { label: 'Bass80', data: bass80_values, borderColor: '#3498db', backgroundColor: 'rgba(52,152,219,0.1)', borderWidth: 2, tension: 0.1, pointRadius: 0 },
                        { label: 'MBass', data: mbass80_values, borderColor: '#2ecc71', backgroundColor: 'rgba(46,204,113,0.1)', borderWidth: 2, tension: 0.1, pointRadius: 0 }
                    ]
                },
                options: {
                    responsive: true, maintainAspectRatio: false,
                    plugins: {
                        title: { display: true, text: `Nucleus-Nucleus Potential`, font: { size: 16 } },
                        legend: { display: true, position: 'top' },
                        tooltip: { mode: 'index', intersect: false, callbacks: { label: (ctx) => `${ctx.dataset.label}: ${ctx.parsed.y.toFixed(2)} MeV` } }
                    },
                    scales: {
                        x: { title: { display: true, text: ' Distance R (fm)' }, grid: { color: 'rgba(0,0,0,0.1)' }, ticks: { maxTicksLimit: 15, autoSkip: true } },
                        y: { title: { display: true, text: ' V (MeV)' }, grid: { color: 'rgba(0,0,0,0.1)' } }
                    },
                    interaction: { intersect: false, mode: 'nearest' }
                }
            });
        }

        // 主逻辑
        document.getElementById('calculateBtn').addEventListener('click', function() {
            const A1 = parseFloat(document.getElementById('a1').value);
            const Z1 = parseFloat(document.getElementById('z1').value);
            const A2 = parseFloat(document.getElementById('a2').value);
            const Z2 = parseFloat(document.getElementById('z2').value);

            if (isNaN(A1) || isNaN(Z1) || isNaN(A2) || isNaN(Z2)) {
                alert('Please enter valid numbers for A1, Z1, A2, Z2.');
                return;
            }

            // 计算势垒
            const results = calculateBarrierHeights(A1, Z1, A2, Z2);

            // 显示反应体系
            const elem1 = getElementSymbol(Z1);
            const elem2 = getElementSymbol(Z2);
            document.getElementById('reactionText').textContent = `${A1}${elem1} + ${A2}${elem2}`;

            // 更新势垒表格
            document.getElementById('bw91-vb').textContent = results.bw91.vb.toFixed(2);
            document.getElementById('bw91-rb').textContent = results.bw91.rb.toFixed(2);
            document.getElementById('bass80-vb').textContent = results.bass80.vb.toFixed(2);
            document.getElementById('bass80-rb').textContent = results.bass80.rb.toFixed(2);
            document.getElementById('mws-vb').textContent = results.mws.vb.toFixed(2);
            document.getElementById('mws-rb').textContent = results.mws.rb.toFixed(2);
            document.getElementById('mbass80-vb').textContent = results.mbass80.vb.toFixed(2);
            document.getElementById('mbass80-rb').textContent = results.mbass80.rb.toFixed(2);

            // ---------- 计算 MWS 的 d0 (用于 Bcap) ----------
            const r0 = 1.2;
            const A1_13 = Math.pow(A1, 1/3);
            const A2_13 = Math.pow(A2, 1/3);
            const d0_mws = r0 * (A1_13 + A2_13) ;

            // 计算 Bcap (BW91 和 MWS)
            const bw91_bcap = computeBcap(BW91, A1, Z1, A2, Z2, d0_mws);
            const mws_bcap = computeBcap(VMWS, A1, Z1, A2, Z2, d0_mws);

            document.getElementById('bw91-bcap').textContent = bw91_bcap.toFixed(2);
            document.getElementById('mws-bcap').textContent = mws_bcap.toFixed(2);

            // 显示结果区域
            document.getElementById('results').style.display = 'block';

            // 绘制势能曲线
            plotPotentialCurves(A1, Z1, A2, Z2, results, `${A1}${elem1} + ${A2}${elem2}`);
        });

        // 初始加载时自动计算一次 (默认16O+208Pb)
        window.addEventListener('DOMContentLoaded', () => {
            document.getElementById('calculateBtn').click();
        });
    })();
</script>
</body>
</html>