FEKO_ISAR二维成像_1

Hsupering

已于 2024-09-18 11:49:07 修改

阅读量229

点赞数 10

分类专栏： FEKO 文章标签： lua

于 2024-09-11 23:47:51 首次发布

本文链接：https://blog.csdn.net/Hsupering/article/details/142142386

版权

FEKO 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

示例：实现对简易物体ISAR二维成像

0000、

由于项目原因接触到了FEKO，正在尝试学习FEKO成像

在校空闲时间会一直更新 : )

一、生成变量

1、将LUA脚本文件拖入CADFEKO，点击运行

2、输入两个维度的分辨率、范围和中心频率

3、脚本会自动生成对应的频率和角度参数

代码附录A

--部分LUA代码 
function getResults()
    local xr = xRange.Value
    local yr = yRange.Value
    local dx = xResolution.Value
    local dy = yResolution.Value
    local f0 = frequencyBase.Value
    local f0n = f0/cf.Const.c0
    local nx = 2*math.ceil(xr/(2*dx))+1
    local ny = 2*math.ceil(yr/(2*dy))+1
    local fxrp = 0.5/dx
    local fyrp = 0.5/dy
    local fxr = fxrp*(nx-1)/nx
    local fyr = fyrp*(ny-1)/ny
    local f1n = f0n - fxr/2
    local phi1r = math.deg(2 * math.atan2(fyr, 2*f1n))
    local f2n = math.sqrt((fxr+f1n)^2 + (fyr/2)^2)
    local f2 = f2n*cf.Const.c0
    local f1 = f1n*cf.Const.c0

    vf1 = ({pcall(function () return nil,project.Variables["f1"] end)})[3] or project.Variables:Add("f1",0)
    vf2 = ({pcall(function () return nil,project.Variables["f2"] end)})[3] or project.Variables:Add("f2",0)
    vnf = ({pcall(function () return nil,project.Variables["nf"] end)})[3] or project.Variables:Add("nf",0)
    vdf = ({pcall(function () return nil,project.Variables["df"] end)})[3] or project.Variables:Add("df",0)
    vRange = ({pcall(function () return nil,project.Variables["range"] end)})[3] or project.Variables:Add("range",0)
    vt1 = ({pcall(function () return nil,project.Variables["phi1"] end)})[3] or project.Variables:Add("phi1",0)
    vt2 = ({pcall(function () return nil,project.Variables["phi2"] end)})[3] or project.Variables:Add("phi2",0)
    vdt = ({pcall(function () return nil,project.Variables["dphi"] end)})[3] or project.Variables:Add("dphi",0)
    vf = ({pcall(function () return nil,project.Variables["f0"] end)})[3] or project.Variables:Add("f0",0)
    vrx = ({pcall(function () return nil,project.Variables["xRange"] end)})[3] or project.Variables:Add("xRange",0)
    vry = ({pcall(function () return nil,project.Variables["yRange"] end)})[3] or project.Variables:Add("yRange",0)
    vdx = ({pcall(function () return nil,project.Variables["xres"] end)})[3] or project.Variables:Add("xres",0)
    vdy = ({pcall(function () return nil,project.Variables["yres"] end)})[3] or project.Variables:Add("yres",0)
end

二、二维ISAR成像

1、设置远场、网格划分，运行FEKO Solver

2、进度条跑完后打开POSTFEKO

3、将LUA脚本拖入POSTFEKO，运行

4、输入相应参数，点击OK

5、输出ISAR二维成像

代码附录B

--部分LUA代码 
local function runForm()
    form:Reject()
    assert(ffSelector.Value ~= nil or ffSelector.Value ~= "", "No farfield result selected")
    ff = ffSelector.Value:GetDataSet()
    assert(ff.Quantities[3].Name == "RCSFactor", "Far Field does not contain RCS Data")
    rForm = pf.Form.New("Angle selection",pf.Enums.FormLayoutEnum.Grid)

    angIndex = 3
    if ff.Axes[angIndex].Count < 2 then
        pf.Form.Critical("Error","Multiple phi angles required for the plane wave source. The script requires variation in the phi direction.")
    end
    if ff.Axes["Frequency"].Count < 2 then
        pf.Form.Critical("Error","The script requires multiple frequency points.")
    end
    local lang0sel = {}
    for i = 2, #ff.Axes[angIndex] - 1 do
        lang0sel[i-1] = ff.Axes[angIndex].Values[i]
    end
    ang0sel = pf.FormComboBox.New("", lang0sel)
    ang0sel.Index = math.ceil(#lang0sel/2)
    local dphi = ff.Axes[angIndex][2] - ff.Axes[angIndex][1]

    rForm:Add(pf.FormLabel.New("Viewing angle:"), 1, 1)
    rForm:Add(ang0sel, 1, 2)   
    
    local function getInputandCalc()
        rForm:Reject()
        freqStart = ff.Axes[1][1]
        freqEnd = ff.Axes[1][#ff.Axes[1]]
        freqSamples = #ff.Axes[1]

        ang1Index = ang0sel.Index - angrsel.Index + 1
        ang2Index = ang0sel.Index + angrsel.Index + 1
        ang1 = ff.Axes[angIndex][ang1Index]
        ang2 = ff.Axes[angIndex][ang2Index]
        ang0 = tonumber(ang0sel.Value)
        angStart = ang1 - ang0
        angEnd = ang2 - ang0
        angSamples = ang2Index - ang1Index + 1
        function num2name(num, idp)
            if math.abs(num) < tonumber("1e-"..idp) then num = 0. end
            return string.gsub(string.gsub(string.format("%." .. (idp or 0) .. "f", num),"[.]","p"),"-","m")
        end
        angrdeg = math.abs(ang2-ang1)
        angrName = num2name(math.abs(ang2-ang1),2)

        f0 = (freqStart + freqEnd)/2
        fd = (freqEnd - freqStart)/(freqSamples-1)
        fRange = fd * freqSamples

        ad = (angEnd - angStart)/(angSamples-1)
        angRange = ad * angSamples

        nFreq = math.floor(freqSamples * upsampleFactor.Value + 0.5)
        nAng = math.floor(angSamples * upsampleFactor.Value + 0.5)

        rx = pf.Const.c0/(2*fRange)
        ry = pf.Const.c0/(4*f0*math.sin(math.rad(angRange)/2))

        nFreqf = math.floor(nFreq/2)
        nAngf = math.floor(nAng/2)

        if osSelector.Value == "Theta" then Sampler = ThetaSampler
        elseif osSelector.Value == "Phi" then Sampler = PhiSampler
        elseif osSelector.Value == "RHC" then Sampler = RHCSampler
        elseif osSelector.Value == "LHC" then Sampler = LHCSampler end

        function SamplerXY(x, y)
            local freq = math.sqrt(x^2 + y^2)
            local ang = math.deg(math.atan2(y, x))

            local freqSample = clamp((freq-freqStart)/fRange*freqSamples, 0, freqSamples-1) + 1
            local angSample = clamp((ang-angStart)/angRange*angSamples, 0, angSamples-1) + 1 + ang1Index - 1
            
            local freqSampleL = math.floor(freqSample)
            local freqSampleU = math.ceil(freqSample)
            local freqSampleD = freqSample - freqSampleL
            local angSampleL = math.floor(angSample)
            local angSampleU = math.ceil(angSample)
            local angSampleD = angSample - angSampleL
            
            local ll = freqSampleD * angSampleD
            local ul = (1 - freqSampleD) * angSampleD
            local uu = (1 - freqSampleD) * (1 - angSampleD)
            local lu = freqSampleD * (1 - angSampleD)
            
            return uu*Sampler(freqSampleL, angSampleL) + ul*Sampler(freqSampleL, angSampleU) +
                   lu*Sampler(freqSampleU, angSampleL) + ll*Sampler(freqSampleU, angSampleU)
        end

        m = pf.ComplexMatrix.New(nFreq, nAng, 0*j)

        fxStart = freqStart
        fyStart = math.tan(math.rad(angStart))*fxStart
        fyEnd = math.tan(math.rad(angEnd))*fxStart
        fxEnd = math.sqrt(freqEnd^2-fyStart^2)
        ysd = (fyEnd - fyStart)/(angSamples-1)
        xsd = (fxEnd - fxStart)/(freqSamples-1)

        if iCheckBox.Checked then
            for f=1,freqSamples do
                for a=1,angSamples do
                    m[f][a] = SamplerXY(fxStart + (f-1)*xsd, fyStart + (a-1)*ysd)
                end
            end
        else
            local f1s = 1
            local f2s = #ff.Axes[1]
            local a1s = ang1Index
            local a2s = ang2Index
            for f=f1s,f2s do
                for a=a1s, a2s do
                    m[f - f1s + 1][a - a1s + 1] = Sampler(f, a)
                end
            end
        end

        if     wfSelector.Value == "Dirichlet Window" then 
            wf = DirichletWindow
        elseif wfSelector.Value == "Hamming Window" then
            wf = HammingWindow
        elseif wfSelector.Value == "Bartlett Window" then
            wf = BartlettWindow
        elseif wfSelector.Value == "Blackman Window" then
            wf = BlackmanWindow 
        end

        m = ApplyWindowMatrix(m, wf, freqSamples, angSamples)
        mOut =FFT2D(m)/(freqSamples*angSamples)
        m = RotateRight(mOut, nFreqf, nAngf)

        if iCheckBox.Checked then
            dfx = xsd/pf.Const.c0
            dfy = ysd/pf.Const.c0
            dr1c = 1/(2 * nFreq*dfx)
            dr2c = 1/(2 * nAng*dfy)
            xRange = 2 * nFreqf * dr1c
            yRange = 2 * nAngf * dr2c
            xStart = -xRange/2
            yStart = -yRange/2
            xEnd = dr1c * (nFreq-1) + xStart
            yEnd = dr2c * (nAng-1)  + yStart
        else
            drx = freqSamples/ nFreq * rx
            dry = angSamples/ nAng * ry
            yRange = 2 * nAngf * dry
            xRange = 2* nFreqf * drx
            xStart = -xRange/2
            yStart = -yRange/2
            xEnd = drx * (nFreq-1) + xStart
            yEnd = dry * (nAng-1)  + yStart
        end

        ds = pf.DataSet.New()
        ds.Axes:Add("frequency", "Hz", f0)
        ds.Axes:Add("X", "m", xStart, xEnd, nFreq)
        ds.Axes:Add("Y", "m", yStart, yEnd, nAng)
        ds.Axes:Add("Z","m",zOffset.Value)
        ds.Quantities:Add("ISAR_Field", "complex", "V")
        ds.Quantities:Add("ISAR_RCS", "scalar", "m^2")
        phi0 = math.rad(ang0)

        rotMat = pf.Matrix.New(3,3,0)
        cphi = math.cos(phi0)
        sphi = math.sin(phi0)
        cphi1 = 1 - cphi

        u1 = ff.MetaData.UVector[1]
        u2 = ff.MetaData.UVector[2]
        u3 = ff.MetaData.UVector[3]
        uVec = pf.Matrix.New(3,1,0)
        uVec[1][1] = u1
        uVec[2][1] = u2
        uVec[3][1] = u3

        v1 = ff.MetaData.VVector[1]
        v2 = ff.MetaData.VVector[2]
        v3 = ff.MetaData.VVector[3]
        vVec = pf.Matrix.New(3,1,0)
        vVec[1][1] = v1
        vVec[2][1] = v2
        vVec[3][1] = v3

        n1 = u2*v3 - u3*v2
        n2 = u3*v1 - u1*v3
        n3 = u1*v2 - u2*v1

        phiv = ang0 + math.deg(math.atan2(n2, n1))
        thtv = math.deg(math.atan2(math.sqrt(n1^2 + n2^2), n3))

        uvec1 = rotMat*uVec
        vvec1 = rotMat*vVec

        ds.MetaData.UVector = pf.Point(uvec1[1][1], uvec1[2][1], uvec1[3][1])
        ds.MetaData.VVector = pf.Point(vvec1[1][1], vvec1[2][1], vvec1[3][1])
        ds.MetaData.Origin = pf.Point(0,0,0)

        for f=1,nFreq do
            for p = 1,nAng do
                ds[1][f][p][1].ISAR_Field = m[f][p]
                ds[1][f][p][1].ISAR_RCS = 4 * pf.Const.pi * Complex.Abs(m[f][p])^2 / magnitude^2
            end
        end

        resdat = ds:StoreData("Custom")
        phi0deg = math.deg(phi0)
        phi0name = num2name(phi0deg,2)
        DatasetName = "ISAR_" .. osSelector.Value .. "_" .. wfSelector.Value:sub(1,wfSelector.Value:find(" ")-1) .. (iCheckBox.Checked and "_P2C" or "").. "_s".. upsampleFactor.Value .. "_v" ..  phi0name.. "_ar" .. angrName
        resdat.Label = nextName(DatasetName, app.StoredData)

        if display3DView.Checked then
            if (ffSelector.Value.Configuration ~= nil) then
                view = app.Views:Add(ffSelector.Value.Configuration)
                plot = view.Plots:Add(resdat)
                view.WindowTitle = resdat.Label
                plot.Visualisation.Opacity = 60
                plot.Quantity.ValuesScaledToDB = true
                plot.Quantity.Type = "ISAR_RCS"
                view.Format.PhiDirection = phiv
                view.Format.ThetaDirection = thtv
                view:ZoomToExtents()
            else
                print("No configuration associated with the RCS data. Please create 3D view and add stored ISAR data manually")
            end
        end

        function MinMaxSampler(i, j)
            return ds[1][i][j][1].ISAR_RCS
        end
        function MinMaxMap(i, j)
            return clamp(i, 1, nFreq), clamp(j, 1, nAng)
        end

        if calcMax.Checked then
            print("\nMaxima")
            maxs = localExtrema(MinMaxSampler, 1, nFreq, 1, nAng, nil, 0, 1000000, MinMaxMap)
            for i=1,maxNum.Value do
                if maxs[i] then
                    print(string.format("Max %s: Value = %s, Xdr = %s, Ycr = %s", i, maxs[i].value, ds.Axes[2][maxs[i].region[1].x], ds.Axes[3][maxs[i].region[1].y]))
                end
            end
        end

        filename = resdat.Label .. ".csv"
        if writeCSV.Checked then
            assert(calcMax.Checked, "Did not calcualte Maximums. Nothing to write to CSV!")
            
            local file = assert(io.open( filename, "w" ), "Could not create the file for writing. Ensure that you have write access.")
            file:write("ISAR")
            file:write("\nModel, " .. ffSelector.Value.Configuration.Model.Label .. "." .. ffSelector.Value.Configuration.Label .. "." .. ffSelector.Value.Label)
            file:write("\nComponent, " .. osSelector.Value)
            file:write("\nWindowing function, " .. wfSelector.Value)
            file:write("\nInterpolated, " .. ( iCheckBox.Checked and "Yes" or "No"))
            file:write("\nView angle [deg], " .. phi0deg)
            file:write("\nAngle range [deg], " .. angrdeg)
            file:write("\n\nMaximums\nIndex, Value, X, Y\n")
            for i=1,maxNum.Value do
                if maxs[i] then
                    file:write(string.format("%s, %s, %s, %s\n", i, maxs[i].value, ds.Axes[2][maxs[i].region[1].x], ds.Axes[3][maxs[i].region[1].y]))
                end
            end
            file:close()
        end
    end
    local function viewingRange()
        rForm:Reject()
        local nphi = math.min(ang0sel.Index,#ff.Axes[angIndex]-ang0sel.Index-1)
        local lRanges = {}
        for i=1,nphi do
            lRanges[i] = i*2*dphi
        end
        angrsel = pf.FormComboBox.New("", lRanges)
        angrsel.Index = #lRanges
        rForm = pf.Form.New("Angle selection",pf.Enums.FormLayoutEnum.Grid)
        rForm:Add(pf.FormLabel.New("Angle range:"), 1, 1)
        rForm:Add(angrsel, 1, 2)
        rForm.Buttons.OK:SetCallBack(getInputandCalc)
        rForm:Run()
    end
    rForm.Buttons.OK:SetCallBack(viewingRange)
    rForm:Run()
end