修改主函数部分使得输入的半径是一个范围，然后输出对应的各个半径消光截面、散射截面、吸收截面、单散射反照率

为了实现主函数中输入的半径是一个范围，并输出对应各个半径的消光截面、散射截面、吸收截面和单散射反照率，可以按照以下步骤进行修改：

1. 修改输入文件格式，使其支持半径范围。
2. 在主函数中读取半径范围，并在循环中调用 `SPHERE_SCAT` 子程序。
3. 将结果写入输出文件。

以下是修改后的代码示例：

### 输入文件格式
假设输入文件 `INPUT.txt` 格式如下：

&input
lamada = 0.5
refre = 1.33
refim = 0.0008
radius_min = 0.5
radius_max = 2.0
radius_step = 0.1
/

### 修改后的主函数

program main
    implicit none
    real, parameter :: pi = 3.141592653589793
    INTEGER, PARAMETER :: NPNMAX = 200000
    real*8 :: lamada, refre, refim, radius_min, radius_max, radius_step, radius
    real*8 :: X, QEXT, QSCA, QBACK, QAB, SSA
    integer :: NANG, NMAX, i
    COMPLEX :: AN_MIE(NPNMAX), BN_MIE(NPNMAX)
    namelist /input/ lamada, refre, refim, radius_min, radius_max, radius_step

    open(unit=13, file='INPUT.txt')
    read(13, nml=input)

    NANG = 91
    NMAX = 10

    open(unit=14, file='output.txt')

    do while (radius_min <= radius_max)
        radius = radius_min
        X = 2 * pi * radius / lamada

        call SPHERE_SCAT(NMAX, lamada, radius, refre, refim, NANG, QEXT, QSCA, QBACK, AN_MIE, BN_MIE)

        QAB = QEXT - QSCA
        SSA = QSCA / QEXT

        write(14, "(F10.6, 4F15.6)") radius, QEXT * pi * radius**2, QSCA * pi * radius**2, QAB, SSA

        radius_min = radius_min + radius_step
    end do

    close(14)
    close(13)
end program

! ********************* 此程序为主要程序 ************************!
!** 函数名：SPHERE_SCAT(NMAX,LAM,RAD,REFRE,REFIM,NANG,QEXT,QSCA,&
! QBACK,AN_MIE,BN_MIE)!
!** 功能：计算：单个粒子的散射矩阵、散射截面、消光截面和吸收截面以及后向散射截面!
!** 参数说明： (输入参数)!
!** NMAX : 球Bessel函数的最大展开阶数!
!** LAM ：入射光波长!
!** RAD : 粒子半径!
!** REFREL: 相对复折射率!
!** NANG : 0-90度散射角区间的角度个数!
!** 参数说明：（输出参数)!
!** S1,S2 : 散射振幅（复向量)!
!** QEXT ：消光效率因子!
!** QSCA : 散射效率因子!
!** QBACK : 后向散射效率因子!
!** AN_MIE: Mie散射系数!
!** BN_MIE: Mie散射系数!

SUBROUTINE SPHERE_SCAT(NMAX, LAM, RAD, REFRE, REFIM, NANG, QEXT, QSCA, QAB, AN_MIE, BN_MIE)
    implicit none
    real, parameter :: pi = 3.141592653589793
    INTEGER, PARAMETER :: NPNMAX = 200000
    real*8, intent(in) :: RAD, LAM, REFRE, REFIM
    integer, intent(in) :: NANG
    INTEGER, INTENT(INOUT) :: NMAX
    COMPLEX :: REFREL, S1(200), S2(200)
    real*8 :: S11(2*NANG-1), S12(2*NANG-1)
    real*8 :: S33(2*NANG-1), S34(2*NANG-1)
    real*8, intent(out) :: QEXT, QSCA, QAB, QBACK
    real*8 :: ANG(2*NANG-1), K0
    COMPLEX :: AN_MIE(NPNMAX), BN_MIE(NPNMAX)
    real :: REFMED
    real :: DANG
    integer :: J, AJ, NAN

    REFMED = 1.0
    K0 = 2 * pi / LAM
    REFREL = CMPLX(REFRE, REFIM) / REFMED
    DANG = pi / FLOAT(NANG - 1) / 2

    CALL LORENZ_MIE(NMAX, LAM, RAD, REFREL, NANG, S1, S2, QEXT, QSCA, QBACK, AN_MIE, BN_MIE)

    QEXT = QEXT * pi * RAD**2
    QSCA = QSCA * pi * RAD**2
    QAB = QEXT - QSCA

    return
end subroutine SPHERE_SCAT

!** 函数名：LORENZ_MIE (NMAX,LAM,RAD,REFREL,NANG,S1,S2,QEXT,QSCA,!
! & QBACK,AN_MIE,BN_MIE)!
!** 功 能：计算球形粒子的散射振幅矩阵，消光、散射及后向散射效率因子，an及bn!
!** 参数说明： (输入参数)!
!** NMAX : 球Bessel函数的最大展开阶数!
!** LAM ：入射光波长!
!** RAD : 粒子半径!
!** REFREL: 相对复折射率!
!** NANG : 0-90度散射角区间的角度个数!
!** 参数说明：（输出参数)!
!** S1,S2 : 散射振幅（复向量)!
!** QEXT ：消光效率因子!
!** QSCA : 散射效率因子!
!** QBACK : 后向散射效率因子!
!** AN_MIE: Mie散射系数!
!** BN_MIE: Mie散射系数!

SUBROUTINE LORENZ_MIE(NMAX, LAM, RAD, REFREL, NANG, S1, S2, QEXT, QSCA, QBACK, AN_MIE, BN_MIE)
    implicit none
    INTEGER, PARAMETER :: NPNMAX = 200000
    REAL*8, INTENT(IN) :: LAM, RAD
    COMPLEX, INTENT(IN) :: REFREL
    INTEGER, INTENT(IN) :: NANG
    INTEGER, INTENT(INOUT) :: NMAX
    COMPLEX, INTENT(OUT) :: S1(200), S2(200)
    REAL*8, INTENT(OUT) :: QEXT, QSCA, QBACK
    COMPLEX, INTENT(OUT) :: AN_MIE(NPNMAX), BN_MIE(NPNMAX)
    REAL*8 :: AMU(100), THETA(100), P1(100), TAU(100), P10(100), P11(100)
    REAL*8 :: X, PI
    COMPLEX :: D(3000000), Y, X1, X10, X11, AN, BN
    DOUBLE PRECISION PS10, PS11, PS1, DN, DX
    INTEGER :: J, NN, N, RN, FN, JJ

    PI = DACOS(-1D0)
    X = 2 * PI * RAD / LAM
    DX = X
    Y = X * REFREL

    XSTOP = X + 4. * X**0.3333 + 2.0
    IF (NMAX >= XSTOP) THEN
        NSTOP = NMAX
    ELSE
        NSTOP = XSTOP
        NMAX = XSTOP
    END IF

    YMOD = CABS(Y)
    NMX = MAX(XSTOP, YMOD) + 15
    DANG = 3.141592653589793 / 2 / FLOAT(NANG - 1)

    DO J = 1, NANG
        THETA(J) = (FLOAT(J) - 1.) * DANG
        AMU(J) = COS(THETA(J))
    END DO

    D(NMX) = CMPLX(0.0, 0.0)
    NN = NMX - 1
    DO N = 1, NN
        RN = NMX - N + 1
        D(NMX - N) = (RN / Y) - (1. / (D(NMX - N + 1) + RN / Y))
    END DO

    DO J = 1, NANG
        P10(J) = 0.0
        P11(J) = 1.0
    END DO

    NN = 2 * NANG - 1
    DO J = 1, NN
        S1(J) = CMPLX(0.0, 0.0)
        S2(J) = CMPLX(0.0, 0.0)
    END DO

    PS10 = DCOS(DX)
    PS11 = DSIN(DX)
    CH10 = -SIN(X)
    CH11 = COS(X)
    APS10 = PS10
    APS11 = PS11
    X10 = CMPLX(APS10, -CH10)
    X11 = CMPLX(APS11, -CH11)
    QSCA = 0.0

    N = 1
    DO WHILE (N - 1 < NSTOP)
        DN = N
        RN = N
        FN = (2. * RN + 1.) / (RN * (RN + 1.))
        PS1 = (2. * DN - 1.) * PS11 / DX - PS10
        APS1 = PS1
        CH1 = (2. * RN - 1.) * CH11 / X - CH10
        X1 = CMPLX(APS1, -CH1)
        AN = (D(N) / REFREL + RN / X) * APS1 - APS11
        AN = AN / ((D(N) / REFREL + RN / X) * X1 - X11)
        BN = (REFREL * D(N) + RN / X) * APS1 - APS11
        BN = BN / ((REFREL * D(N) + RN / X) * X1 - X11)
        AN_MIE(N) = AN
        BN_MIE(N) = BN
        QSCA = QSCA + (2. * RN + 1.) * (CABS(AN)**2 + CABS(BN)**2)

        DO J = 1, NANG
            JJ = 2 * NANG - J
            P1(J) = P11(J)
            TAU(J) = RN * AMU(J) * P1(J) - (RN + 1.) * P10(J)
            P = (-1.)**(N - 1)
            S1(J) = S1(J) + FN * (AN * P1(J) + BN * TAU(J))
            T = (-1.)**N
            S2(J) = S2(J) + FN * (AN * TAU(J) + BN * P1(J))
            IF (J == JJ) THEN
                CYCLE
            END IF
            S1(JJ) = S1(JJ) + FN * (AN * P1(J) * P + BN * TAU(J) * T)
            S2(JJ) = S2(JJ) + FN * (AN * TAU(J) * T + BN * P1(J) * P)
        END DO

        PS10 = PS11
        PS11 = PS1
        APS11 = PS11
        CH10 = CH11
        CH11 = CH1
        X11 = CMPLX(APS11, -CH11)
        N = N + 1
        RN = N
        DO J = 1, NANG
            P11(J) = ((2. * RN - 1.)
            P10(J) = P1(J)
        END DO
    END DO

    QSCA = (2. / (X * X)) * QSCA
    QEXT = (4. / (X * X)) * REAL(S1(1))
    QBACK = (4. / (X * X)) * CABS(S1(2 * NANG - 1))**2

    return
end subroutine LORENZ_MIE

### 输出文件格式
输出文件 `output.txt` 将包含每一半径下的消光截面、散射截面、吸收截面和单散射反照率，格式如下：

radius   extinction_cross_section   scattering_cross_section   absorption_cross_section   single_scattering_albedo
0.500000       0.000000           0.000000           0.000000           0.000000
0.600000       0.000000           0.000000           0.000000           0.000000
...

通过上述修改，程序将能够处理半径范围内的多个值，并输出每个半径对应的物理量。