import numpy as np
import math
import sys
import os

# --- 定数の定義 ---
MAXDATA_SIZE = 101
JIKU_SIZE    = 4
DATA1_SIZE   = 7
DATA_SIZE    = 21
DATA2_SIZE   = 31

PI_CONST = math.pi

# --- グローバル変数 (NumPy配列として定義、1ベースインデックス対応のためサイズ+1) ---
InFile = ""
OutFile = ""
in_stream = None
out_stream = None
SilentMode = 0

x = np.zeros(DATA_SIZE + 1, dtype=np.float64)
sx = np.zeros(DATA_SIZE + 1, dtype=np.float64)

ai = np.zeros(JIKU_SIZE + 1, dtype=np.float64)
dai = np.zeros(JIKU_SIZE + 1, dtype=np.float64)
ao = np.zeros(JIKU_SIZE + 1, dtype=np.float64)
dao = np.zeros(JIKU_SIZE + 1, dtype=np.float64)
cao = np.zeros(JIKU_SIZE + 1, dtype=np.float64)
dcao = np.zeros(JIKU_SIZE + 1, dtype=np.float64)
agi = np.zeros(JIKU_SIZE + 1, dtype=np.float64)
dagi = np.zeros(JIKU_SIZE + 1, dtype=np.float64)
ago = np.zeros(JIKU_SIZE + 1, dtype=np.float64)
dago = np.zeros(JIKU_SIZE + 1, dtype=np.float64)

cai = np.zeros(DATA1_SIZE + 1, dtype=np.float64)
dcai = np.zeros(DATA1_SIZE + 1, dtype=np.float64)
sai = np.zeros(DATA1_SIZE + 1, dtype=np.float64)
dsai = np.zeros(DATA1_SIZE + 1, dtype=np.float64)
sao = np.zeros(DATA1_SIZE + 1, dtype=np.float64)
dsao = np.zeros(DATA1_SIZE + 1, dtype=np.float64)

vi = 0.0
dvi = 0.0
vo = 0.0
dvo = 0.0

p = np.zeros(DATA_SIZE + 1, dtype=np.float32)
sigp = np.zeros(DATA_SIZE + 1, dtype=np.float32)

bo = np.zeros(MAXDATA_SIZE + 1, dtype=np.float32)
wgt = np.zeros(MAXDATA_SIZE + 1, dtype=np.float32)
xlsq = np.zeros((DATA_SIZE + 1, MAXDATA_SIZE + 1), dtype=np.float32)
bc = np.zeros(MAXDATA_SIZE + 1, dtype=np.float32)
bd = np.zeros(MAXDATA_SIZE + 1, dtype=np.float32)
sigb = np.zeros(MAXDATA_SIZE + 1, dtype=np.float32)
bw = np.zeros(MAXDATA_SIZE + 1, dtype=np.float32)
fh = np.zeros((JIKU_SIZE + 1, MAXDATA_SIZE + 1), dtype=np.float32)

iinum = np.zeros(MAXDATA_SIZE + 1, dtype=int)
ih1 = np.zeros(MAXDATA_SIZE + 1, dtype=int)
ih2 = np.zeros(MAXDATA_SIZE + 1, dtype=int)
ih3 = np.zeros(MAXDATA_SIZE + 1, dtype=int)

am = np.zeros((DATA2_SIZE + 1, DATA2_SIZE + 1), dtype=np.float32)
aam = np.zeros((DATA2_SIZE + 1, DATA2_SIZE + 1), dtype=np.float32)

# --- 関数定義 ---

# ## 関数: volume
def volume():
    term = 1.0 - x[4] * x[4] - x[5] * x[5] - x[6] * x[6] + 2.0 * x[4] * x[5] * x[6]
    if term < 0:
        return 0.0
    return (x[1] * x[2] * x[3] * math.sqrt(term))

# ## 関数: sinal
def sinal():
    term = 1.0 - x[1] * x[1]
    if term < 0:
        return 0.0
    return (math.sqrt(term))

# ## 関数: acosd
def acosd():
    arg = x[1]
    if arg < -1.0: arg = -1.0
    if arg > 1.0: arg = 1.0

    if arg == 0.0:
        ac = 90.0
    else:
        ac = math.acos(arg) * 57.29578
    return (ac)

# ## 関数: cosaif
def cosaif():
    if x[2] == 0.0 or x[3] == 0.0:
        return 0.0
    return (x[1] / (x[2] * x[3]))

# ## 関数: trigcf
def trigcf():
    denominator = (1.0 - x[1] * x[1])
    if denominator == 0.0:
        return 0.0
    return (x[1] * (x[1] - 1.0) / denominator)

# ## 関数: trigaf
def trigaf():
    denominator_val = (1.0 - 3.0 * x[2] * x[2] + 2.0 * x[2] * x[2] * x[2])
    if denominator_val <= 0.0 or x[1] == 0.0:
        return 0.0
    numerator_val = (1.0 - x[2] * x[2])
    if numerator_val < 0.0:
        return 0.0
    return (math.sqrt(numerator_val / denominator_val) / x[1])

# ## 関数: trigvf
def trigvf():
    val = 1.0 - 3.0 * x[2] * x[2] + 2.0 * x[2] * x[2] * x[2]
    if val < 0:
        return 0.0
    return (x[1] * x[1] * x[1] * math.sqrt(val))

# ## 関数: sigmaf
def sigmaf(y_val, n, f_func):
    global x, sx

    sigy_sq = 0.0
    
    current_y = f_func() 
    y_val[0] = current_y 

    for i in range(1, n + 1):
        xt = x[i] 
        
        if abs(sx[i]) < 1e-10: 
            continue 
        
        x[i] = 0.0001 * sx[i] + x[i]
        
        ff = (f_func() - current_y) * 10000.0 
        sigy_sq += ff * ff 
        x[i] = xt 
    
    return math.sqrt(sigy_sq)

# ## 関数: inv30s (行列の逆行列計算)
def inv30s(a_matrix, b_matrix, n):
    try:
        temp_a = a_matrix[1:n+1, 1:n+1].astype(np.float64) 
        
        inv_temp_a = np.linalg.inv(temp_a)
        
        b_matrix[1:n+1, 1:n+1] = inv_temp_a.astype(np.float32)

    except np.linalg.LinAlgError:
        out_stream.write("No Solution !! (Singular matrix)\n")
        sys.exit(-1)

# ## 関数: lsq (最小二乗法)
# no: number of observations
# num_params: number of parameters (renamed from 'np')
def lsq(no, num_params, wl): # Renamed 'np' to 'num_params'
    global p, sigp, bo, wgt, xlsq, bc, bd, sigb, bw, iinum, ih1, ih2, ih3, am, aam

    vt = np.zeros(DATA_SIZE + 1, dtype=np.float32)
    oe = np.zeros(MAXDATA_SIZE + 1, dtype=np.float32)

    aaa = float(no - num_params) # Use num_params

    while True:
        if aaa < 0:
            out_stream.write("Number of observation is smaller than that of parameters.\n")
            out_stream.write("The problem can not be solved\n")
            return -1
        elif aaa == 0:
            out_stream.write("Number of observations is equal to that of parameters.\n")
            out_stream.write("Parameters are obtained without least square calculation.\n")
            
            am[1:num_params+1, 1:num_params+1] = xlsq[1:num_params+1, 1:num_params+1].copy() # Use num_params

            inv30s(am, aam, num_params) # Use num_params
            for i in range(1, num_params + 1): # Use num_params
                p[i] = 0.0
                for k in range(1, num_params + 1): # Use num_params
                    p[i] += aam[k][i] * bo[k]
            return no
        else:
            for i in range(1, num_params + 1): # Use num_params
                vt[i] = 0.0
                for k in range(1, num_params + 1): # Use num_params
                    am[i][k] = 0.0
                    for j in range(1, no + 1):
                        am[i][k] += xlsq[i][j] * wgt[j] * xlsq[k][j]
                for j in range(1, no + 1):
                    vt[i] += xlsq[i][j] * wgt[j] * bo[j]
                if am[i][i] <= 0:
                    out_stream.write("Normal equation contains zero or negative diagonalelement.\n")
                    out_stream.write("The problem can not be solved.\n")
                    return -1
            
            inv30s(am, aam, num_params) # Use num_params
            
            for i in range(1, num_params + 1): # Use num_params
                p[i] = 0.0
                for k in range(1, num_params + 1): # Use num_params
                    p[i] += aam[k][i] * vt[k]
            
            sig = 0.0
            for j in range(1, no + 1):
                bc[j] = 0.0
                for i in range(1, num_params + 1): # Use num_params
                    bc[j] += xlsq[i][j] * p[i]
                sig += wgt[j] * (bo[j] - bc[j]) * (bo[j] - bc[j])
            
            sqrsig = math.sqrt(sig / aaa)
            ot = 0.0
            for j in range(1, no + 1):
                if wgt[j] <= 0.0:
                    sigb[j] = 0.0
                else:
                    sigb[j] = sqrsig / math.sqrt(wgt[j])
                
                bd[j] = bo[j] - bc[j]
                
                if abs(bd[j]) - abs(3.0 * sigb[j]) < 0 :
                    oe[j] = 0
                else:
                    oe[j] = 1.0
                    ot += 1.0
                    
                    arg_asin = wl * math.sqrt(max(0.0, bc[j])) / 2.0
                    if arg_asin < -1.0 or arg_asin > 1.0:
                        btc = 0.0
                    else:
                        btc = 2.0 * math.asin(arg_asin)
                    
                    btc /= 0.0174533
                    dif = bw[j] - btc
                    out_stream.write("This observation is eliminated.\n")
                    out_stream.write(f"No. = {j:2d}  bo = {bw[j]:10.3f}  bc = {btc:10.3f} dif(B) = {dif:10.3f}\n")
            
            if ot != 0:
                k = 1
                for j in range(1, no + 1):
                    if oe[j] == 0:
                        bo[k] = bo[j]
                        wgt[k] = wgt[j]
                        iinum[k] = iinum[j]
                        ih1[k] = ih1[j]
                        ih2[k] = ih2[j]
                        ih3[k] = ih3[j]
                        bw[k] = bw[j]
                        for i in range(1, num_params + 1): # Use num_params
                            xlsq[i][k] = xlsq[i][j]
                        k += 1
                no = k - 1
                aaa = float(no - num_params) # Use num_params
            else:
                for i in range(1, num_params + 1): # Use num_params
                    sigp[i] = math.sqrt(aam[i][i]) * sqrsig
                return no

# ## 関数: recip
def recip():
    global x, sx, ai, dai, ao, dao, cao, dcao, agi, dagi, ago, dago, cai, dcai, sai, dsai, sao, dsao, vi, dvi, vo, dvo

    n = 1
    y_temp = np.zeros(1, dtype=np.float64) 

    for j in range(1, 3 + 1):
        x[1] = cai[j]
        sx[1] = dcai[j]
        dsai[j] = sigmaf(y_temp, n, sinal)
        sai[j] = y_temp[0] 
        cai[j + 3] = cai[j]
        dcai[j + 3] = dcai[j]
        sai[j + 3] = sai[j]
        dsai[j + 3] = dsai[j]
    
    m = 3
    for j in range(1, 3 + 1):
        for k in range(1, 3 + 1):
            jk = j + k
            x[k] = cai[jk - 1] 
            sx[k] = dcai[jk - 1]
        dcao[j] = sigmaf(y_temp, m, rco) 
        cao[j] = y_temp[0] 
        
        x[1] = cao[j]
        sx[1] = dcao[j]
        dsao[j] = sigmaf(y_temp, n, sinal)
        sao[j] = y_temp[0] 
        sao[j + 3] = sao[j]
        dsao[j + 3] = dsao[j]
    
    n = 4
    for j in range(1, 3 + 1):
        x[1] = ai[j]
        x[2] = cai[j]
        x[3] = cai[j + 1]
        x[4] = cai[j + 2]
        sx[1] = dai[j]
        sx[2] = dcai[j]
        sx[3] = dcai[j + 1]
        sx[4] = dcai[j + 2]
        dao[j] = sigmaf(y_temp, n, rax) 
        ao[j] = y_temp[0] 

    for j in range(1, 3 + 1):
        x[j] = ai[j]
        x[j + 3] = cai[j] 
        sx[j] = dai[j]
        sx[j + 3] = dcai[j]
    
    n = 6
    dvi = sigmaf(y_temp, n, volume) 
    vi = y_temp[0] 
    
    if vi == 0.0:
        vo = 0.0
        dvo = 0.0
    else:
        vo = 1.0 / vi
        dvo = dvi * vo * vo
    
    n = 1
    for j in range(1, 3 + 1):
        x[1] = cai[j]
        sx[1] = dcai[j]
        dagi[j] = sigmaf(y_temp, n, acosd) 
        agi[j] = y_temp[0] 
        
        x[1] = cao[j]
        sx[1] = dcao[j]
        dago[j] = sigmaf(y_temp, n, acosd) 
        ago[j] = y_temp[0] 


# ## メイン関数
def main(argv):
    global InFile, OutFile, in_stream, out_stream, SilentMode
    global x, sx, ai, dai, ao, dao, cao, dcao, agi, dagi, ago, dago
    global cai, dcai, sai, dsai, sao, dsao
    global vi, dvi, vo, dvo
    global p, sigp, bo, wgt, xlsq, bc, bd, sigb, bw, fh
    global iinum, ih1, ih2, ih3
    global am, aam

    le = np.zeros(DATA1_SIZE + 1, dtype=int)
    ih = np.zeros(JIKU_SIZE + 1, dtype=int)
    xh = np.zeros(DATA1_SIZE + 1, dtype=np.float32)

    title_buffer = ""

    dobs = 0.0
    dcal = 0.0
    thecc = 0.0
    thecl = 0.0
    theta = 0.0

    # Corrected: Added num_params back to the unpack list and initialized with 0
    i, ie, io, ip, iw, j, jj, job, jw, k, ls, n, no, num, num_params = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
    abc, def_val, ghi, pqr, xyz, bas, wl, wls, wwl, radius, sigg, zo = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

    m, l = 0, 0

    argc = len(argv)
    if argc < 3:
        sys.stderr.write("Usage: K.py InputFile OutputFile [--silent]\n")
        sys.exit(1)
    
    InFile = argv[1]
    OutFile = argv[2]
    if argc >= 4 and "--silent" in argv[3:]:
        SilentMode = 1

    print()
    print(f"InFile: {InFile}")
    print(f"OutFile: {OutFile}")
    
    try:
        in_stream = open(InFile, 'r')
    except FileNotFoundError:
        sys.stderr.write(f"Can not read {InFile}\n")
        sys.exit(2)

    try:
        out_stream = open(OutFile, 'w')
    except IOError:
        sys.stderr.write(f"Can not write to {OutFile}\n")
        in_stream.close()
        sys.exit(2)

    if not SilentMode:
        sys.stderr.write("Least Square Calculation for Lattice Constant\n")
        sys.stderr.write("            modified by T.Kamiya    1989.10\n")
        sys.stderr.write("            modified by T.Kamiya    2005.08.02\n\n")
    out_stream.write("Least Square Calculation for Lattice Constant\n")
    out_stream.write("            modified by T.Kamiya    1989.10\n")
    out_stream.write("            modified by T.Kamiya    2005.08.02\n\n")

    while True:
        line = in_stream.readline()
        if not line:
            out_stream.write("\n\nAbnormal Program END for EOF\n")
            sys.exit(-2)
        title_buffer = line.strip()
        out_stream.write(f"\n {title_buffer}\n")

        try:
            params = in_stream.readline().split()
            ls, le[1], le[2], le[3], le[4], le[5], iw, io, ip = map(int, params)
        except ValueError:
            sys.stderr.write("Error reading main parameters. Invalid format.\n")
            sys.exit(3)

        out_stream.write("      HLS   LE(1)  LE(2)  LE(3)  LE(4)  LE(5)    IW     IO     IP\n")
        out_stream.write(f"{ls:7d}")
        for i in range(1, 5 + 1):
            out_stream.write(f"{le[i]:7d}")
        out_stream.write(f"{iw:7d}{io:7d}{ip:7d}\n")

        num = 0
        jj = 0
        jw = 0

        no = 1
        try:
            wl, radius = map(float, in_stream.readline().split())
        except ValueError:
            sys.stderr.write("Error reading wavelength/radius. Invalid format.\n")
            sys.exit(3)

        out_stream.write(f"Wave Length ={wl:f}         Radius ={radius:f}\n")
        if jw == 0:
            wwl = wl
        jw += 1
        wls = 4.0 / (wl * wl)

        try:
            ih_vals = in_stream.readline().split()
            ih[1], ih[2], ih[3], bo[no] = int(ih_vals[0]), int(ih_vals[1]), int(ih_vals[2]), float(ih_vals[3])
        except ValueError:
            sys.stderr.write("Error reading initial data point. Invalid format.\n")
            sys.exit(3)
        
        # Changed peek() to check for an empty string, which is what readline() returns at EOF.
        # peek() reads the next character without advancing the stream, so it's good for checking if more data exists.
        # But for readline-based parsing, checking the result of readline is more direct for EOF.
        # The loop condition relies on the *next* ih[1] value, which is read at the end of the loop body.
        # The loop must continue as long as `in_stream.good()` and `ih[1]` condition is met.
        # The current structure: read line, then check if it's EOF, then process.
        # The while loop condition `ih[1] < 1000 or ih[1] == 2000` applies to the *last read* ih[1].
        # It means we've already read one set of ih, bo[no] and are checking if we should process it.
        # This loop structure is very common in C programs with `scanf` followed by `while` loop condition.
        # In Python, a do-while like structure or a loop that breaks internally is often used.
        # The `in_stream.peek(1) != ''` is for when the file might end *before* an expected line.
        # However, the subsequent `in_stream.readline()` would return `''` if EOF, and `split()` would then fail.
        # The current error handling for `readline()` returning empty string is likely sufficient.

        while (ih[1] < 1000 or ih[1] == 2000): # Condition for the current ih[1]
            if ih[1] == 2000:
                try:
                    wl, radius = map(float, in_stream.readline().split())
                except ValueError:
                    sys.stderr.write("Error reading new wavelength/radius. Invalid format.\n")
                    sys.exit(3)
                out_stream.write(f"Wave Length ={wl:f}         Radius ={radius:f}\n")
                if jw == 0:
                    wwl = wl
                jw += 1
                wls = 4.0 / (wl * wl)
            else:
                num += 1
                jj += 1
                iinum[jj] = num
                ih1[jj] = ih[1]
                ih2[jj] = ih[2]
                ih3[jj] = ih[3]
                bw[jj] = bo[no]
                
                if io == 1:
                    bo[no] = math.sin(0.01745329 * bo[no])
                elif io == 2:
                    pass
                elif io == 3:
                    try:
                        zo = float(in_stream.readline().strip())
                    except ValueError:
                        sys.stderr.write("Error reading zo for io=3. Invalid format.\n")
                        sys.exit(3)
                    xyz = math.cos(bo[no] / radius)
                    abc = math.atan(zo / radius)
                    def_val = math.cos(abc)
                    ghi = xyz * def_val
                    if ghi < -1.0: ghi = -1.0
                    if ghi > 1.0: ghi = 1.0
                    pqr = math.acos(ghi)
                    bo[no] = math.sin(0.5 * pqr)
                elif io == 4:
                    bo[no] = 0.5 * bo[no]
                    bo[no] = math.sin(0.01745329 * bo[no])
                elif io == 5:
                    if bo[no] == 0.0: bo[no] = 0.0
                    else: bo[no] = wl * 0.5 / bo[no]
                else:
                    abc = bo[no] / radius
                    xyz = 1.0 + abc * abc
                    if xyz <= 0.0: bo[no] = 0.0
                    else: bo[no] = abc / math.sqrt(xyz)
                
                bas = bo[no] * bo[no]

                if iw == 0:
                    try:
                        sigg = float(in_stream.readline().strip())
                    except ValueError:
                        sys.stderr.write("Error reading sigg for iw=0. Invalid format.\n")
                        sys.exit(3)
                    sigg *= 0.25
                    if sigg <= 0: sigg = 0.25
                    if bas * (1.0 - bas) == 0.0: wgt[no] = 0.0
                    else: wgt[no] = sigg / (bas * (1.0 - bas))
                elif iw == 1:
                    try:
                        sigg = float(in_stream.readline().strip())
                    except ValueError:
                        sys.stderr.write("Error reading sigg for iw=1. Invalid format.\n")
                        sys.exit(3)
                    sigg *= 2.0
                    if bas * (1.0 - bas) * sigg * sigg == 0.0: wgt[no] = 0.0
                    else: wgt[no] = 0.25 / (bas * (1.0 - bas) * sigg * sigg)
                elif iw == 2:
                    wgt[no] = 1.0
                elif iw == 3:
                    try:
                        sigg = float(in_stream.readline().strip())
                    except ValueError:
                        sys.stderr.write("Error reading sigg for iw=3. Invalid format.\n")
                        sys.exit(3)
                    if sigg > 0: wgt[no] = sigg
                    else: wgt[no] = 1.0
                elif iw == 4:
                    try:
                        sigg = float(in_stream.readline().strip())
                    except ValueError:
                        sys.stderr.write("Error reading sigg for iw=4. Invalid format.\n")
                        sys.exit(3)
                    if sigg == 0.0: wgt[no] = 0.0
                    else: wgt[no] = 1.0 / (sigg * sigg)
                
                for i in range(1, 3 + 1):
                    aaih_val = float(ih[i])
                    fh[i][no] = aaih_val
                    xh[i] = fh[i][no] * fh[i][no]
                xh[4] = 2.0 * fh[2][no] * fh[3][no]
                xh[5] = 2.0 * fh[3][no] * fh[1][no]
                xh[6] = 2.0 * fh[1][no] * fh[2][no]

                if ls == 1:
                    num_params = 6
                    for i in range(1, 6 + 1): xlsq[i][no] = xh[i]
                elif ls == 2:
                    num_params = 4
                    for i in range(1, 3 + 1): xlsq[i][no] = xh[i]
                    xlsq[4][no] = xh[5]
                elif ls == 3:
                    num_params = 4
                    for i in range(1, 3 + 1): xlsq[i][no] = xh[i]
                    xlsq[4][no] = xh[6]
                elif ls == 4:
                    num_params = 3
                    for i in range(1, 3 + 1): xlsq[i][no] = xh[i]
                elif ls == 5:
                    num_params = 2
                    xlsq[1][no] = xh[1] + xh[2]
                    xlsq[2][no] = xh[3]
                elif ls == 6:
                    num_params = 1
                    xlsq[1][no] = xh[1] + xh[2] + xh[3]
                elif ls == 7:
                    num_params = 2
                    xlsq[1][no] = xh[1] + xh[2] + xh[3]
                    xlsq[2][no] = xh[4] + xh[5] + xh[6]
                elif ls == 8:
                    num_params = 2
                    xlsq[1][no] = xh[1] + xh[2] + 0.5 * xh[6]
                    xlsq[2][no] = xh[3]

                if le[1]: 
                    num_params += 1
                    arg_asin = bo[no]
                    if arg_asin < -1.0: arg_asin = -1.0
                    if arg_asin > 1.0: arg_asin = 1.0
                    
                    val_asin = math.asin(arg_asin)
                    xlsq[num_params][no] = (PI_CONST/2.0 - val_asin) * math.tan(PI_CONST/2.0 - val_asin) 
                if le[2]: num_params += 1; xlsq[num_params][no] = 4.0 * bas * (1.0 - bas)
                if le[3]: 
                    num_params += 1
                    arg_asin = bo[no]
                    if arg_asin < -1.0: arg_asin = -1.0
                    if arg_asin > 1.0: arg_asin = 1.0
                    val_asin = math.asin(arg_asin)
                    if bo[no] == 0.0 or val_asin == 0.0: xlsq[num_params][no] = 0.0
                    else: xlsq[num_params][no] = 4.0 * bas * (1.0 - bas) * (1.0 / bo[no] + 1.0 / val_asin) 
                if le[4]: 
                    num_params += 1
                    arg_asin = bo[no]
                    if arg_asin < -1.0: arg_asin = -1.0
                    if arg_asin > 1.0: arg_asin = 1.0
                    xlsq[num_params][no] = math.sin(2.0 * math.asin(arg_asin)) 
                if le[5]: 
                    num_params += 1
                    arg_asin = bo[no]
                    if arg_asin < -1.0: arg_asin = -1.0
                    if arg_asin > 1.0: arg_asin = 1.0
                    theta = math.asin(arg_asin) 
                    xlsq[num_params][no] = math.sin(2.0 * theta) * math.cos(theta)
                
                bo[no] = bas * wls
                no += 1

            next_line = in_stream.readline()
            if not next_line:
                out_stream.write("\n\nAbnormal Program END for EOF\n")
                sys.exit(-2)
            try:
                ih_vals = next_line.split()
                ih[1], ih[2], ih[3], bo[no] = int(ih_vals[0]), int(ih_vals[1]), int(ih_vals[2]), float(ih_vals[3])
            except ValueError:
                sys.stderr.write("Error reading next data point. Invalid format.\n")
                sys.exit(3)
        
        no -= 1 
        no = lsq(no, num_params, wl)

        out_stream.write("            P(I)                 SIGP(I)\n")
        for i in range(1, num_params + 1):
            out_stream.write(f"          {p[i]:10.7f}           {sigp[i]:10.7f}\n")
        
        out_stream.write("          H     K     L    Dobs.    Dcal.     BO      BC        D(B)      d(2Q)\n")
        for j in range(1, no + 1):
            sqrt_bo_j = math.sqrt(max(0.0, bo[j]))
            dobs = 1.0 / sqrt_bo_j if sqrt_bo_j != 0 else float('inf')

            sqrt_bc_j = math.sqrt(max(0.0, bc[j]))
            dcal = 1.0 / sqrt_bc_j if sqrt_bc_j != 0 else float('inf')

            thecc = wwl / (2.0 * dcal) if dcal != 0 else float('inf')
            
            arg_asin_thecc = thecc
            if arg_asin_thecc < -1.0: arg_asin_thecc = -1.0
            if arg_asin_thecc > 1.0: arg_asin_thecc = 1.0
            
            thecl = 2.0 * math.asin(arg_asin_thecc) * 57.29578
            
            out_stream.write(f" ({iinum[j]:2d}){ih1[j]:5d}{ih2[j]:5d}{ih3[j]:5d}{dobs:8.4f}{dcal:8.4f}{bw[j]:8.3f}{thecl:8.3f}  {bd[j]:9.5f}  {(bw[j]-thecl):9.5f}\n")

        for i in range(1, 3 + 1):
            cai[i] = 0.0; dcai[i] = 0.0
            cao[i] = 0.0; dcao[i] = 0.0
            agi[i] = 90.0; dagi[i] = 0.0
            ago[i] = 90.0; dago[i] = 0.0

        aai = np.zeros(DATA1_SIZE + 1, dtype=np.float32)
        daai = np.zeros(DATA1_SIZE + 1, dtype=np.float32)

        if ls == 1:
            y_temp_sig = np.zeros(1, dtype=np.float64)

            for i in range(1, 3 + 1):
                ai[i] = math.sqrt(float(p[i])) if p[i] >= 0 else 0.0
                if ai[i] == 0.0: dai[i] = 0.0
                else: dai[i] = sigp[i] / (2.0 * ai[i])
                
                aai[i] = float(ai[i])
                aai[i + 3] = float(ai[i])
                daai[i] = float(dai[i])
                daai[i + 3] = float(dai[i])
            
            n = 3
            for i in range(1, 3 + 1):
                x[1] = p[i + 3]
                sx[1] = sigp[i + 3]
                for j in range(1, 2 + 1):
                    k = i + j
                    x[j + 1] = aai[k]
                    sx[j + 1] = daai[k]
                
                dcai[i] = sigmaf(y_temp_sig, n, cosaif)
                cai[i] = y_temp_sig[0]
            recip()
        elif ls == 2:
            p[5] = p[4]; sigp[5] = sigp[4]
            p[4] = 0.0; sigp[4] = 0.0
            p[6] = 0.0; sigp[6] = 0.0
            
            y_temp_sig = np.zeros(1, dtype=np.float64)

            for i in range(1, 3 + 1):
                ai[i] = math.sqrt(float(p[i])) if p[i] >= 0 else 0.0
                if ai[i] == 0.0: dai[i] = 0.0
                else: dai[i] = sigp[i] / (2.0 * ai[i])
            
            for i in range(1, 3 + 1):
                aai[i] = float(ai[i])
                aai[i + 3] = float(ai[i])
                daai[i] = float(dai[i])
                daai[i + 3] = float(dai[i])
            n = 3
            for i in range(1, 3 + 1):
                x[1] = p[i + 3]; sx[1] = sigp[i + 3]
                for j in range(1, 2 + 1):
                    k = i + j
                    x[j + 1] = aai[k]; sx[j + 1] = daai[k]
                dcai[i] = sigmaf(y_temp_sig, n, cosaif)
                cai[i] = y_temp_sig[0]
            recip()
        elif ls == 3:
            p[6] = p[4]; sigp[6] = sigp[4]
            p[4] = 0.0; sigp[4] = 0.0
            p[5] = 0.0; sigp[5] = 0.0
            
            y_temp_sig = np.zeros(1, dtype=np.float64)

            for i in range(1, 3 + 1):
                ai[i] = math.sqrt(float(p[i])) if p[i] >= 0 else 0.0
                if ai[i] == 0.0: dai[i] = 0.0
                else: dai[i] = sigp[i] / (2.0 * ai[i])
            
            for i in range(1, 3 + 1):
                aai[i] = float(ai[i])
                aai[i + 3] = float(ai[i])
                daai[i] = float(dai[i])
                daai[i + 3] = float(dai[i])
            n = 3
            for i in range(1, 3 + 1):
                x[1] = p[i + 3]; sx[1] = sigp[i + 3]
                for j in range(1, 2 + 1):
                    k = i + j
                    x[j + 1] = aai[k]; sx[j + 1] = daai[k]
                dcai[i] = sigmaf(y_temp_sig, n, cosaif)
                cai[i] = y_temp_sig[0]
            recip()
        elif ls == 4:
            for i in range(1, 3 + 1):
                ai[i] = math.sqrt(float(p[i])) if p[i] >= 0 else 0.0
                if ai[i] == 0.0: dai[i] = 0.0
                else: dai[i] = sigp[i] / (2.0 * ai[i])
                if ai[i] == 0.0: ao[i] = 0.0
                else: ao[i] = 1.0 / ai[i]
                dao[i] = dai[i] * ao[i] * ao[i]
            
            if any(val == 0.0 for val in [ai[1], ai[2], ai[3]]): vi = 0.0
            else: vi = ai[1] * ai[2] * ai[3]
            
            div_term1 = (dai[1] / ai[1])**2 if ai[1] != 0 else 0.0
            div_term2 = (dai[2] / ai[2])**2 if ai[2] != 0 else 0.0
            div_term3 = (dai[3] / ai[3])**2 if ai[3] != 0 else 0.0

            term_dvi = div_term1 + div_term2 + div_term3
            if term_dvi < 0: dvi = 0.0
            else: dvi = vi * math.sqrt(term_dvi)
            
            if vi == 0.0: vo = 0.0; dvo = 0.0
            else: vo = 1.0 / vi; dvo = dvi * vo * vo
        elif ls == 5:
            ai[3] = math.sqrt(float(p[2])) if p[2] >= 0 else 0.0
            if ai[3] == 0.0: dai[3] = 0.0
            else: dai[3] = 0.5 * sigp[2] / ai[3]
            ai[1] = math.sqrt(float(p[1])) if p[1] >= 0 else 0.0
            if ai[1] == 0.0: dai[1] = 0.0
            else: dai[1] = 0.5 * sigp[1] / ai[1]
            ai[2] = ai[1]
            dai[2] = dai[1]
            for j in range(1, 3 + 1):
                if ai[j] == 0.0: ao[j] = 0.0
                else: ao[j] = 1.0 / ai[j]
                dao[j] = dai[j] * ao[j] * ao[j]
            
            if any(val == 0.0 for val in [ai[3], ai[1]]): vi = 0.0
            else: vi = ai[3] * ai[1] * ai[1]
            
            div_term1 = (dai[1] / ai[1])**2 if ai[1] != 0 else 0.0
            div_term3 = (dai[3] / ai[3])**2 if ai[3] != 0 else 0.0
            term_dvi = 2.0 * div_term1 + div_term3
            
            if term_dvi < 0: dvi = 0.0
            else: dvi = vi * math.sqrt(term_dvi)
            
            if vi == 0.0: vo = 0.0; dvo = 0.0
            else: vo = 1.0 / vi; dvo = dvi * vo * vo
        elif ls == 6:
            ai[1] = math.sqrt(float(p[1])) if p[1] >= 0 else 0.0
            if ai[1] == 0.0: dai[1] = 0.0
            else: dai[1] = 0.5 * sigp[1] / ai[1]
            if ai[1] == 0.0: ao[1] = 0.0
            else: ao[1] = 1.0 / ai[1]
            dao[1] = dai[1] * ao[1] * ao[1]
            for j in range(2, 3 + 1):
                ai[j] = ai[1]; dai[j] = dai[1]
                ao[j] = ao[1]; dao[j] = dao[1]
            
            vi = ai[1] * ai[1] * ai[1]
            dvi = 3.0 * ai[1] * ai[1] * dai[1]
            
            if vi == 0.0: vo = 0.0; dvo = 0.0
            else: vo = 1.0 / vi; dvo = dvi * vo * vo
        elif ls == 7:
            y_temp_sig = np.zeros(1, dtype=np.float64)
            for j in range(1, 3 + 1):
                ai[j] = math.sqrt(float(p[1])) if p[1] >= 0 else 0.0
                if p[1] == 0.0: cai[j] = 0.0
                else: cai[j] = p[2] / p[1]
                if ai[1] == 0.0: dai[j] = 0.0
                else: dai[j] = 0.5 * sigp[1] / ai[1]
                
                div_term1 = (sigp[1] / p[1])**2 if p[1] != 0 else 0.0
                div_term2 = (sigp[2] / p[2])**2 if p[2] != 0 else 0.0
                term_dcai = div_term1 + div_term2
                if p[1] == 0.0 or p[2] == 0.0 or term_dcai < 0: dcai[j] = 0.0
                else: dcai[j] = cai[j] * math.sqrt(term_dcai)
            
            x[1] = ai[1]; sx[1] = dai[1]
            x[2] = cai[1]; sx[2] = dcai[1]
            n = 2
            dao[1] = sigmaf(y_temp_sig, n, trigaf)
            ao[1] = y_temp_sig[0]

            dvi = sigmaf(y_temp_sig, n, trigvf)
            vi = y_temp_sig[0]
            
            x[1] = cai[1]; sx[1] = dcai[1]
            dcao[1] = sigmaf(y_temp_sig, n, trigcf)
            cao[1] = y_temp_sig[0]

            dagi[1] = sigmaf(y_temp_sig, n, acosd)
            agi[1] = y_temp_sig[0]
            
            x[1] = cao[1]; sx[1] = dcao[1]
            dago[1] = sigmaf(y_temp_sig, n, acosd)
            ago[1] = y_temp_sig[0]
            
            for j in range(2, 3 + 1):
                cao[j] = cao[1]; ago[j] = ago[1]; agi[j] = agi[1]
                dagi[j] = dagi[1]; dcao[j] = dcao[1]; dago[j] = dago[1]
                dao[j] = dao[1]; ao[j] = ao[1]
            
            if vi == 0.0: vo = 0.0; dvo = 0.0
            else: vo = 1.0 / vi; dvo = dvi * vo * vo
        elif ls == 8:
            ai[1] = math.sqrt(float(p[1])) if p[1] >= 0 else 0.0
            ai[2] = ai[1]
            ai[3] = math.sqrt(float(p[2])) if p[2] >= 0 else 0.0
            
            if ai[1] == 0.0: dai[1] = 0.0
            else: dai[1] = 0.5 * sigp[1] / ai[1]
            dai[2] = dai[1]
            if ai[3] == 0.0: dai[3] = 0.0
            else: dai[3] = 0.5 * sigp[2] / ai[3]
            
            if any(val == 0.0 for val in [ai[1], ai[2], ai[3]]): vi = 0.0
            else: vi = 0.8660254 * ai[1] * ai[2] * ai[3]
            
            div_term1 = (dai[1] / ai[1])**2 if ai[1] != 0 else 0.0
            div_term3 = (dai[3] / ai[3])**2 if ai[3] != 0 else 0.0
            term_dvi = 4.0 * div_term1 + div_term3
            if term_dvi < 0: dvi = 0.0
            else: dvi = vi * math.sqrt(term_dvi)
            
            if ai[1] == 0.0: ao[1] = 0.0
            else: ao[1] = 1.1547005 / ai[1]
            ao[2] = ao[1]
            if ai[3] == 0.0: ao[3] = 0.0
            else: ao[3] = 1.0 / ai[3]
            
            if ai[1] == 0.0: dao[1] = 0.0
            else: dao[1] = dai[1] * ao[1] / ai[1]
            dao[2] = dao[1]
            if ai[3] == 0.0: dao[3] = 0.0
            else: dao[3] = dai[3] * ao[3] / ai[3]
            
            cai[3] = 0.5
            cao[3] = -0.5
            agi[3] = 60.0
            ago[3] = 120.0
            
            if vi == 0.0: vo = 0.0; dvo = 0.0
            else: vo = 1.0 / vi; dvo = dvi * vo * vo
        
        out_stream.write("Reciprocal cell constant\n")
        out_stream.write("      a* b* c*\n")
        for j in range(1, 3 + 1):
            out_stream.write(f"{ai[j]:10.7f}({dai[j]:10.7f})    ")
        out_stream.write("\n    alpha* beta* gamma*\n")
        for j in range(1, 3 + 1):
            out_stream.write(f"{agi[j]:10.6f}({dagi[j]:10.6f})    ")
        out_stream.write("\n    cos a* cos b* cos c*\n")
        for j in range(1, 3 + 1):
            out_stream.write(f"{cai[j]:10.6f}({dcai[j]:10.6f})    ")
        out_stream.write(f"\n    V*={vi:.10g}({dvi:.10g})\n")

        out_stream.write("\nDirect cell constant\n")
        out_stream.write("        a                         b                         c\n")
        for j in range(1, 3 + 1):
            out_stream.write(f"{ao[j]:8.5f}({dao[j]:8.5f})        ")
        out_stream.write("\n      alpha                     beta                      gamma\n")
        for j in range(1, 3 + 1):
            out_stream.write(f"{ago[j]:8.4f}({dago[j]:8.4f})        ")
        out_stream.write("\n      cos a                     cos b                     cos c\n")
        for j in range(1, 3 + 1):
            out_stream.write(f"{cao[j]:8.4f}({dcao[j]:8.4f})        ")
        out_stream.write(f"\n    V={vo:.10g}({dvo:.10g})\n")

        next_line = in_stream.readline()
        if not next_line:
            out_stream.write("\n\nAbnormal Program END for EOF\n")
            sys.exit(-2)
        try:
            job = int(next_line.strip())
        except ValueError:
            sys.stderr.write("Error reading job code. Invalid format.\n")
            sys.exit(3)
        
        if job == 0:
            break

    if not SilentMode: out_stream.write("\nLattice Constant Calculation was finished!\n")
    print("\nLattice Constant Calculation was finished!\n")

    in_stream.close()
    out_stream.close()

    print(f"LSQ results are saved to [{OutFile}]")
    input("\nPress ENTER to terminate>>\n")
    
    sys.exit(0)

if __name__ == "__main__":
    main(sys.argv)