Optimizing controls deflection on a flying wing...

[joyo44]

Hey Peter and followers,

I'm trying to find the optimized (CD = min) controls deflection for a swept flying wing and a given design point( best glide ratio or climb ratio), the wing has 3 controls per side.

I'm in doubt about the results from ASB (neuralfoil), specifically when I compare them to AVL or AVL+Xfoil.

For best climb ratio I got the following results from ASB:

Interesting are the low deflections for a dCL of 0.195, the shift of nearly 3% in static margin and a dAlpha of 2.6°.
CL = 0.609 is the point for best climb found by ASB.

If I use AVL with these deflections the results look like this:

which gives dCL = 0.03, SM shift = 0 % and dAlpha = 0.33°.

With AVL+Xfoil:

which gives dCL = 0.031, SM shift = 0.1 % and dAlpha = 0.57°. The best climb ratio for this calculation is around CL = 0.64 and not
within reach with these deflections.

I don't know where these discrepancies come from, is it a shortcoming of ASB/Neuralfoil or is my setup of the script wrong.

My code is as follows:

#//////////////////////////////////////////////////////////////////////////////
# Optimization: Flying Wing with Flaps
#//////////////////////////////////////////////////////////////////////////////
#
# ASB_FlyingWingOptFlap.py
#
#
# JY, 01.2024    ASB 4.1.7
#//////////////////////////////////////////////////////////////////////////////
#
#------------------------------------------------------------------------------
import aerosandbox as asb
import aerosandbox.numpy as np
import aerosandbox.tools.units as u
import aerosandbox.tools.pretty_plots as p
from aerosandbox.tools.string_formatting import eng_string

import pandas as pd
import copy
import sys
import os
import json
import datetime
import scipy
import matplotlib.pyplot as plt
import shutil
import argparse

from Utils import *
from Dynamic_Stability import *

#------------------------------------------------------------------------------
IpoptTol = 0.001
g = 9.81
rad = 180.0/np.pi

# Dict to collect airfoil data
AirfoilData = {"Name": "Undef",
"Thickness": 0.0,
"Camber": 0.0
}
AirfoilDataList = []

# Get Date & Time  
now = datetime.datetime.now()

WorkingDir = os.getcwd()
RootDir = os.path.abspath('..')
PostFix = ""

###############################################################################
##### Section: Input  #########################################################
###############################################################################

Project         = "Opt Flap Deflection"
AeroMethod      = "AeroBuildup" # AeroBuildup (| LLT | NLLT --> not used)
Neuralfoil_size = "large"       # medium | large | xxxlarge

#------------------------------------------------------------------------------
# Enironment
Altitude    = 0.0 # m
ISAday      = 0.0 # dT to ISA

#------------------------------------------------------------------------------
# OperatingPoint
TrimCM  = False
TrimGLZ = False
TrimSTZ = True
TrimCL  = False
CLtrim  = 0.5
TrimAlpha = False
Alphatrim = 5.0

if TrimCM:
    PostFix = "TrimCM"       
elif TrimGLZ:
    PostFix = "TrimGLZ"      
elif TrimSTZ:
    PostFix = "TrimSTZ"       
elif TrimCL:
    PostFix = "TrimCL"       
elif TrimAlpha:
    PostFix = "TrimAlpha" 

#------------------------------------------------------------------------------
# Model
Span = 2.8        # m  
Model_Mass = 2.0  # kg 
Model_CoG  =  [0.29764, 0.0, 0.0]

#------------------------------------------------------------------------------
# Wing Segments             
Wing_Spanwise_panel_resolution  = [10]
Wing_Iwing     = 0.0
Wing_LEpos     = [[    0.0,       0.0, 0.0],
			      [0.20573, 0.4666667, 0.0],
			      [0.40382, 0.9333333, 0.0],
			      [0.60190, 1.4000000, 0.0]
		         ]
Wing_Chords    = [0.25, 0.23333, 0.21667, 0.2]
Wing_Twists    = [0.0, -0.653, -1.911, -4.22]     
Wing_Airfoils  = ["YO_120-a2","YO_120-a2", "YO_120-a2", "YO_120-a2"]

UseWingControls    = True
UseFixedDeflection = False
if TrimCM:
    UseWingControls = False
    
Wing_control_surface_is_symmetrics = ["True",  "True",  "True"]
Wing_control_names                 = ["flap", "flap1", "elevon"]
Wing_control_surface_deflections   = [   0.0,     0.0,      0.0]      
Wing_control_surface_hinge_points  = [  0.75,    0.75,   0.75,   0.75]

#------------------------------------------------------------------------------
# Airfoil Setup
for i in range(len(Wing_Airfoils)):
    print("i : ", i , "  Airfoil : ", Wing_Airfoils[i])
    af = asb.Airfoil(
        name=Wing_Airfoils[i],
        ).repanel(n_points_per_side=100)  
        
    afW.append(af)
    
#------------------------------------------------------------------------------
opti = asb.Opti()

###############################################################################
##### Section: Wing Build #####################################################
###############################################################################

HalfSpan = Span/2.0
NoOfSections = len(Wing_Chords)
NoOfWingSegments = NoOfSections - 1
XsecsW = []
iSeg = 0     

if UseWingControls:
    #
    # fixed deflection from input
    #
    if UseFixedDeflection:

        for iSect in range(NoOfSections):
            if iSect == NoOfWingSegments: 
                XsecsW.append(asb.WingXSec(
                            xyz_le=Wing_LEpos[iSect],
                            chord=Wing_Chords[iSect],
                            airfoil=afW[iSect],
                            twist=Wing_Twists[iSect],)
                            )    
            else:  
                iSeg = iSect            
                control = asb.ControlSurface(name=Wing_control_names[iSeg],
                                             symmetric=Wing_control_surface_is_symmetrics[iSeg],
                                             trailing_edge=True,
                                             hinge_point=Wing_control_surface_hinge_points[iSect],
                                             deflection=Wing_control_surface_deflections[iSeg],                                             
                                             )   
                XsecsW.append(asb.WingXSec(
                                xyz_le=Wing_LEpos[iSect],
                                chord=Wing_Chords[iSect],
                                airfoil=afW[iSect],
                                twist=Wing_Twists[iSect],
                                control_surfaces=[control],)
                                )     
    else:
        #
        # deflection is an opti_var
        #            
        if TrimCM != True:            
            for iSect in range(NoOfSections): 
                if iSect == NoOfWingSegments: 
                    XsecsW.append(asb.WingXSec(
                                xyz_le=Wing_LEpos[iSect],
                                chord=Wing_Chords[iSect],
                                airfoil=afW[iSect],
                                twist=Wing_Twists[iSect],)
                                )    
                else:  
                    iSeg = iSect 
                
                    if Wing_control_names[iSeg].lower() != 'none':
                        control = asb.ControlSurface(name=Wing_control_names[iSeg],
                                                     symmetric=Wing_control_surface_is_symmetrics[iSeg],
                                                     trailing_edge=True,
                                                     hinge_point=Wing_control_surface_hinge_points[iSect],
                                                     deflection=opti.variable(init_guess=0.0, lower_bound=-10.0, upper_bound=10.0, scale=1)
                                                     )

                        XsecsW.append(asb.WingXSec(
                                        xyz_le=Wing_LEpos[iSect],
                                        chord=Wing_Chords[iSect],
                                        airfoil=afW[iSect],
                                        twist=Wing_Twists[iSect],
                                        control_surfaces=[control],)
                                        )
                    else:
                        XsecsW.append(asb.WingXSec(
                                    xyz_le=Wing_LEpos[iSect],
                                    chord=Wing_Chords[iSect],
                                    airfoil=afW[iSect],
                                    twist=Wing_Twists[iSect],)
                                    ) 
        else:
            for iSect in range(NoOfSections):  
                if iSect == NoOfWingSegments: 
                    XsecsW.append(asb.WingXSec(
                            xyz_le=Wing_LEpos[iSect],
                            chord=Wing_Chords[iSect],
                            airfoil=afW[iSect],
                            twist=Wing_Twists[iSect],)
                            )    
                else:  
                    iSeg = iSect         
                    control = asb.ControlSurface(name=Wing_control_names[iSeg],
                                                 symmetric=Wing_control_surface_is_symmetrics[iSeg],
                                                 trailing_edge=True,
                                                 hinge_point=Wing_control_surface_hinge_points[iSect],
                                                 deflection=0.0,                                             
                                                 )   
                    XsecsW.append(asb.WingXSec(
                                    xyz_le=Wing_LEpos[iSect],
                                    chord=Wing_Chords[iSect],
                                    airfoil=afW[iSect],
                                    twist=Wing_Twists[iSect],
                                    control_surfaces=[control],)
                                    ) 
else:
    for iSect in range(NoOfSections):
            XsecsW.append(asb.WingXSec(
                        xyz_le=Wing_LEpos[iSect],
                        chord=Wing_Chords[iSect],
                        airfoil=afW[iSect],
                        twist=Wing_Twists[iSect],)
                        )

wing = asb.Wing(
    name="Main Wing",
    color='k',
    symmetric=True,
    xsecs=XsecsW,
).subdivide_sections(ratio=Wing_Spanwise_panel_resolution[0], spacing_function=np.cosspace) 
        
for i in range(len(wing.xsecs)):
    print("i = ", i, " xsecs[i] = ", wing.xsecs[i], " xsecs[i].control_surfaces = ", wing.xsecs[i].control_surfaces)
 
###############################################################################
##### Section: Airplane Build #################################################
###############################################################################

airplane = asb.Airplane(
            name="Flying Wing",
            xyz_ref=Model_CoG,
            wings=[wing]
            )

###############################################################################
##### Section: Aerodynamics ###################################################
###############################################################################

Atmosphere = asb.Atmosphere(altitude=Altitude)
Density = Atmosphere.density()
  
velocity=opti.variable(
    init_guess=10.0,
    lower_bound=6.0,
    upper_bound=30.0,
    log_transform=True
)            

if TrimAlpha != True:
    alpha=opti.variable(
        init_guess=5.0,
        lower_bound=-5.0,
        upper_bound=15.0
    )
else:
    alpha = Alphatrim
    
op_point = asb.OperatingPoint(
    atmosphere = Atmosphere,
    velocity=velocity,
    alpha=alpha
    )

#if AeroMethod == "AeroBuildup":
ab = asb.AeroBuildup(
    airplane=airplane,
    op_point=op_point,
    xyz_ref=Model_CoG,
    model_size=Neuralfoil_size
)

aero = ab.run_with_stability_derivatives(
    alpha=True,
    beta=True,
    p=True,
    q=True,
    r=True,
    )

LD = aero['CL']/aero['CD']
CP = aero['CL']**1.5/aero['CD']
Gamma = -np.arctan(aero["CD"]/aero["CL"])*rad
power_loss = aero["D"] * op_point.velocity
sink_rate = power_loss / g / Model_Mass
lift_residual = aero['L'] - Model_Mass * g
static_margin = (aero["x_np"] - Model_CoG[0]) / wing.mean_aerodynamic_chord()
pressure_center = Model_CoG[0] - airplane.c_ref* aero['Cm']/aero['CL']  

if TrimCM:
    print("TrimCM")
    opti.subject_to([
        aero["L"] == g * Model_Mass,          
        aero["Cm"] == 0,      
    ])    

    opti.minimize(lift_residual ** 2)    
elif TrimGLZ:
   print("TrimGLZ")
   opti.subject_to([
        aero["L"] == g * Model_Mass,          
        aero["Cm"] == 0,      
    ])    
   opti.minimize(-LD)  
elif TrimSTZ:
   print("TrimSTZ")
   opti.subject_to([
        aero["L"] == g * Model_Mass,          
        aero["Cm"] == 0,      
    ])    
   opti.minimize(-CP)
elif TrimCL:
   print("TrimCL")
   opti.subject_to([
        aero["L"] == g * Model_Mass,          
        aero["Cm"] == 0,
        aero["CL"] == CLtrim        
    ])    
   opti.minimize(lift_residual ** 2)    
elif TrimAlpha:
   print("TrimAlpha")
   opti.subject_to([
            aero["L"] == g * mass_props_TOGW.mass,     
            op_point.alpha == Alphatrim,
            ])      
   opti.minimize(lift_residual ** 2)   
   
###############################################################################
##### Section: Optimization ###################################################
###############################################################################   
def fmt(x):
    return f"{s(x):.6g}"   
try:
    sol = opti.solve(max_iter=300, behavior_on_failure="return_last", options={'ipopt.tol': IpoptTol})
    s = lambda x: sol.value(x)    
except RuntimeError:
        #------------------------------------------------------------------------------            
        print("#################################### DEBUG ##############################")
        sol = opti.debug
        opti.debug.show_infeasibilities(IpoptTol)
        s = lambda x: sol.value(x)
        print("Velocity        : ", fmt(op_point.velocity))
        print("Gamma           : ", fmt(Gamma))
        print("Sink Rate       : ", fmt(-sink_rate))           
        print("AoA             : ", fmt(op_point.alpha)) 
        print("LD              : ", fmt(LD)) 
        print("CP              : ", fmt(CP))         
        print("CL              : ", fmt(aero["CL"]))
        print("CD              : ", fmt(aero["CD"]))
        print("Cm              : ", fmt(aero["Cm"]))
        print("XNP             : ", fmt(aero["x_np"])) 
        print("XSP             : ", fmt(Model_CoG[0]))
        print("STM             : ", fmt((aero["x_np"] - Model_CoG[0]) / wing.mean_aerodynamic_chord()))

        airplane.draw_three_view(show=False)
        p.show_plot(title=Project, tight_layout=False, savefig="Airplane3view_shaded.png")             
            
        sys.exit() 
        
airplane.substitute_solution(sol)
op_point.substitute_solution(sol)  

aero = {
    k: s(v)
    for k, v in aero.items() if not isinstance(v, list)
}  

#
# Write AVL input
#
wing_lowres = copy.deepcopy(wing)
xsecs_to_keep = np.arange(len(wing.xsecs)) % Wing_Spanwise_panel_resolution[0] == 0
xsecs_to_keep[0] = True
xsecs_to_keep[-1] = True
wing_lowres.xsecs = np.array(wing_lowres.xsecs)[xsecs_to_keep]

airplanAVL = asb.Airplane(
        wings=[wing_lowres]
        ) 

airplanAVL.analysis_specific_options = {
    asb.AVL: dict(
        profile_drag_coefficient=0.011
    )
}

airplanAVL.wings[0].analysis_specific_options = {
    asb.AVL: dict(
        spanwise_resolution=Wing_Spanwise_panel_resolution[0],
        component=1
        )
}

avl_geo = asb.AVL(
    airplane=airplanAVL,
    op_point=op_point,
    xyz_ref=Model_CoG,
    working_directory=WorkingDir,        
)  

avl_geo.write_avl(
                  filepath="airplane.avl"
                  )
###############################################################################
##### Section: Finalize 
###############################################################################

LD = aero["L"] / aero["D"]
Gamma = -np.arctan(aero["CD"]/aero["CL"])*rad
power_loss = aero["D"] * op_point.velocity
sink_rate = power_loss / g / Model_Mass 
    
AOAtrim = op_point.alpha   
CLtrim = aero['CL']
CDtrim = aero['CD']
LDtrim = aero['CL']/aero['CD']
CRtrim = aero['CL']**1.5/aero['CD']
CMtrim = aero['Cm'] 
                       
###############################################################################
##### Section: Polar generation
############################################################################### 
if AeroMethod == "AeroBuildup":
    opti2 = asb.Opti()

    vels = np.linspace(5, 15, 100)

    op_point_sweep = asb.OperatingPoint(
        atmosphere=Atmosphere,
        velocity = vels,
        alpha=opti2.variable(
            init_guess=5 * np.ones_like(vels),
            lower_bound=-5,
            upper_bound=15
        )
    )
    
    aero_sweep = asb.AeroBuildup(
        airplane=airplane,
        op_point=op_point_sweep,
        xyz_ref=Model_CoG,
        model_size=Neuralfoil_size
    ).run()

    lift_residual = aero_sweep['L'] - Model_Mass * g

    opti2.minimize(np.sum(lift_residual ** 2))

    try:
        sol2 = opti2.solve()
    except RuntimeError:
        sol2 = opti2.debug

    s = lambda x: sol2.value(x)
      
    sol2.show_infeasibilities()  

    trimmable_points = (np.abs(sol2.value(lift_residual)) < 1e-3)

    airplane.substitute_solution(sol2)
    op_point.substitute_solution(sol2)  
        
    # Curves
    Alphas = sol2.value(op_point_sweep.alpha)
    CLs    = sol2.value(aero_sweep['CL'])
    CDs    = sol2.value(aero_sweep['CD'])
    CMs    = sol2.value(aero_sweep['Cm'])
    CRs    = sol2.value(aero_sweep['CL']**1.5/aero_sweep['CD'])    
    LDs    = sol2.value(aero_sweep['L'] / aero_sweep['D'])

trimmable_points = CLs > 0.0   

# Plot
airplane.draw_wireframe(
                       color="b",
                       set_axis_visibility=True,
                       show=True,
                       )
                       
fig, ax = plt.subplots(2, 2, figsize=(10, 8))

if TrimCM:
    CalcType = "\nTrim for CM = 0"
elif TrimGLZ:
    CalcType = "\nTrim for best L/D"  
elif TrimSTZ:
    CalcType = "\nTrim for best climb power"
elif TrimCL:
    CalcType = ("\nTrim for CL = %.3f" % (aero['CL'])) 
elif TrimAlpha:
    CalcType = ("\nTrim for Alpha = %.3f°" % (op_point.alpha))        
    
Title1 = Project + "\nAirfoil: " + Wing_Airfoils[0] + CalcType
fig.suptitle(Title1, fontsize=10)

plt.sca(ax[0, 0])
plt.plot(Alphas, CLs)
plt.plot(AOAtrim, CLtrim, color='r', marker='o', markerfacecolor='r', markersize=3)
plt.xlabel("Angle of Attack [deg]")
plt.ylabel("Lift Coefficient $C_L$ [-]")
ax[0, 0].grid(True)

plt.sca(ax[0, 1])
plt.plot(CLs[trimmable_points], LDs[trimmable_points])
plt.plot(CLtrim, LDtrim, color='r', marker='o', markerfacecolor='r', markersize=3)
plt.xlabel("Lift Coefficient $C_L$ [-]")
plt.ylabel("Glide Ratio")
ax[0, 1].grid(True)

plt.sca(ax[1, 0])
plt.plot(CLs[trimmable_points], CRs[trimmable_points])
plt.plot(CLtrim, CRtrim, color='r', marker='o', markerfacecolor='r', markersize=3)
plt.xlabel("Lift Coefficient $C_L$ [-]")
plt.ylabel("Climb Power Ratio")
ax[1, 0].grid(True)

plt.sca(ax[1, 1])
plt.plot(CLs[trimmable_points], CMs[trimmable_points])
plt.plot(CLtrim, CMtrim, color='r', marker='o', markerfacecolor='r', markersize=3)
plt.xlabel("Lift Coefficient $C_L$ [-]")
plt.ylabel("Moment Coefficient $C_M$ [-]")
ax[1, 1].grid(True)

plt.grid(True)

textstr = ""
textstr += (f"XSP = {Model_CoG[0]:.3f} m\n")
textstr += (f"STM = {sol.value(static_margin)*100.0:.1f} %\n")
textstr += (f"V   = {op_point.velocity:.1f} m/s\n")
textstr += (f"Alp = {op_point.alpha:.1f}°\n")    
textstr += (f"CL  = {CLtrim:.2f}\n")   
        
# Get control deflections
LastName = ""
for wi, wing in enumerate(airplane.wings):
    for xi, xsec in enumerate(wing.xsecs):
        for csi, surf in enumerate(xsec.control_surfaces):
                if surf.deflection != 0.0 and surf.name != LastName:
                    textstr += (f"{surf.name} = {surf.deflection:.2f}°\n")
                    LastName = surf.name
                    break
                    

textstr = textstr.rstrip('\n')  

if textstr != "":        
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    fig.text(0.7, 0.3, textstr, transform=ax[0, 0].transAxes, fontsize=8,
        verticalalignment='top', bbox=props)
    
p.show_plot(savefig="AircraftPolar" + "_" + PostFix + ".png",  savefig_transparent=False, show=True)     
    
# Write Polar
f = open("AirplanePolar.out","w+") 
f.write("Trim Point\n")
f.write("==========\n\n")
f.write("Project: %s\n\n" % (Project)) 

f.write("Bref = %.3f m \n" % (airplane.b_ref))
f.write("Cref = %.3f m \n" % (airplane.c_ref))
f.write("Sref = %.3f m² \n\n" % (airplane.s_ref))

f.write("XNP  = %.3f m \n" % (aero["x_np"]))
f.write("XSP  = %.3f m \n" % (Model_CoG[0]))
f.write("STM  = %.3f \n\n" % (sol.value(static_margin))) 

f.write("AoA = %.2f\n" % (AOAtrim))
f.write("CL  = %.3f\n" % (CLtrim))
f.write("CD  = %.5f\n" % (CDtrim)) 
f.write("CM  = %.5f\n" % (CMtrim))
f.write("CM0 = %.5f\n" % (sol.value(static_margin)*CLtrim))     
f.write("L/D = %.2f\n" % (LDtrim))   
f.write("STZ = %.2f\n\n" % (CRtrim))     

# Get control deflections
textstr = ""
LastName = ""
for wi, wing in enumerate(airplane.wings):
    for xi, xsec in enumerate(wing.xsecs):
        for csi, surf in enumerate(xsec.control_surfaces):
                if surf.deflection != 0.0 and surf.name != LastName:
                    textstr += (f"{surf.name} = {surf.deflection:.2f}°\n")
                    LastName = surf.name
                    break
if textstr != "":      
    textstr = textstr.rstrip('\n')  

f.write("%s \n\n" % (textstr))   
    
f.write("  i     Alpha        CL        CD        CM       GLZ       STZ\n")
f.write("---------------------------------------------------------------\n")    

for i in range(len(Alphas)):
    f.write("%3d%10.2f%10.3f%10.4f%10.4f%10.2f%10.2f\n" % (i, Alphas[i], CLs[i], CDs[i], CMs[i], LDs[i], CRs[i]))
              
f.close()        
                      

and the airfoil:

YO_120-a2
    1.000000      0.6798118E-10
   0.9916942      0.3500178E-04
   0.9806790      0.1273933E-03
   0.9682153      0.3342214E-03
   0.9545113      0.6806938E-03
   0.9399032      0.1174611E-02
   0.9247018      0.1817735E-02
   0.9091382      0.2607375E-02
   0.8933557      0.3539338E-02
   0.8774318      0.4608118E-02
   0.8614053      0.5807756E-02
   0.8452917      0.7131421E-02
   0.8290964      0.8571635E-02
   0.8128190      0.1012017E-01
   0.7964575      0.1176792E-01
   0.7800094      0.1350529E-01
   0.7634737      0.1532207E-01
   0.7468500      0.1720761E-01
   0.7301405      0.1915105E-01
   0.7133485      0.2114159E-01
   0.6964810      0.2316779E-01
   0.6795511      0.2521807E-01
   0.6625826      0.2728012E-01
   0.6456165      0.2933967E-01
   0.6286821      0.3138518E-01
   0.6117904      0.3340761E-01
   0.5949439      0.3539999E-01
   0.5781409      0.3735644E-01
   0.5613793      0.3927171E-01
   0.5446558      0.4114127E-01
   0.5279683      0.4296070E-01
   0.5113147      0.4472578E-01
   0.4946934      0.4643278E-01
   0.4781040      0.4807775E-01
   0.4615471      0.4965648E-01
   0.4450234      0.5116520E-01
   0.4285357      0.5259942E-01
   0.4120862      0.5395455E-01
   0.3956798      0.5522576E-01
   0.3793216      0.5640745E-01
   0.3630184      0.5749345E-01
   0.3467772      0.5847696E-01
   0.3306053      0.5935028E-01
   0.3145086      0.6010469E-01
   0.2984885      0.6073125E-01
   0.2825449      0.6122300E-01
   0.2666833      0.6157323E-01
   0.2509133      0.6177482E-01
   0.2352492      0.6181987E-01
   0.2197083      0.6169971E-01
   0.2043147      0.6140527E-01
   0.1890976      0.6092679E-01
   0.1740952      0.6025458E-01
   0.1593580      0.5937953E-01
   0.1449506      0.5829388E-01
   0.1309564      0.5699331E-01
   0.1174789      0.5547920E-01
   0.1046401      0.5376143E-01
   0.9257054E-01  0.5186057E-01
   0.8139122E-01  0.4980814E-01
   0.7119183E-01  0.4764386E-01
   0.6201245E-01  0.4540968E-01
   0.5383922E-01  0.4314351E-01
   0.4661293E-01  0.4087429E-01
   0.4024679E-01  0.3862131E-01
   0.3464307E-01  0.3639433E-01
   0.2970554E-01  0.3419654E-01
   0.2534723E-01  0.3202643E-01
   0.2149316E-01  0.2987978E-01
   0.1808065E-01  0.2775129E-01
   0.1505871E-01  0.2563520E-01
   0.1238668E-01  0.2352607E-01
   0.1003246E-01  0.2141967E-01
   0.7971321E-02  0.1931314E-01
   0.6183914E-02  0.1720549E-01
   0.4654965E-02  0.1509853E-01
   0.3371661E-02  0.1299635E-01
   0.2321719E-02  0.1090619E-01
   0.1492387E-02  0.8837201E-02
   0.8682410E-03  0.6800073E-02
   0.4316667E-03  0.4806004E-02
   0.1560951E-03  0.2865784E-02
   0.2229662E-04  0.9869924E-03
   0.2806187E-04 -0.8360795E-03
   0.4784696E-03 -0.2618314E-02
   0.1483701E-02 -0.4226822E-02
   0.2886974E-02 -0.5626729E-02
   0.4569133E-02 -0.6864648E-02
   0.6481300E-02 -0.7986768E-02
   0.8609820E-02 -0.9027252E-02
   0.1096056E-01 -0.1001066E-01
   0.1355266E-01 -0.1095525E-01
   0.1641634E-01 -0.1187541E-01
   0.1959345E-01 -0.1278277E-01
   0.2313917E-01 -0.1368743E-01
   0.2712448E-01 -0.1459870E-01
   0.3163971E-01 -0.1552507E-01
   0.3679975E-01 -0.1647471E-01
   0.4275303E-01 -0.1745569E-01
   0.4967382E-01 -0.1847348E-01
   0.5774342E-01 -0.1952764E-01
   0.6712960E-01 -0.2061037E-01
   0.7793558E-01 -0.2170403E-01
   0.9013800E-01 -0.2278078E-01
   0.1035736     -0.2380884E-01
   0.1179883     -0.2476153E-01
   0.1331155     -0.2562287E-01
   0.1487337     -0.2638684E-01
   0.1646831     -0.2705470E-01
   0.1808575     -0.2763134E-01
   0.1971877     -0.2812304E-01
   0.2136301     -0.2853602E-01
   0.2301557     -0.2887612E-01
   0.2467450     -0.2914855E-01
   0.2633846     -0.2935770E-01
   0.2800646     -0.2950751E-01
   0.2967778     -0.2960130E-01
   0.3135184     -0.2964195E-01
   0.3302814     -0.2963186E-01
   0.3470642     -0.2957305E-01
   0.3638629     -0.2946750E-01
   0.3806760     -0.2931676E-01
   0.3975008     -0.2912230E-01
   0.4143360     -0.2888517E-01
   0.4311805     -0.2860667E-01
   0.4480326     -0.2828769E-01
   0.4648918     -0.2792910E-01
   0.4817574     -0.2753177E-01
   0.4986288     -0.2709638E-01
   0.5155050     -0.2662380E-01
   0.5323859     -0.2611452E-01
   0.5492716     -0.2556909E-01
   0.5661612     -0.2498842E-01
   0.5830549     -0.2437273E-01
   0.5999528     -0.2372291E-01
   0.6168545     -0.2303946E-01
   0.6337599     -0.2232289E-01
   0.6506691     -0.2157377E-01
   0.6675826     -0.2079274E-01
   0.6844997     -0.1998044E-01
   0.7014208     -0.1913744E-01
   0.7183457     -0.1826432E-01
   0.7352746     -0.1736173E-01
   0.7522081     -0.1643032E-01
   0.7691450     -0.1547083E-01
   0.7860859     -0.1448382E-01
   0.8030296     -0.1347009E-01
   0.8199755     -0.1243038E-01
   0.8369198     -0.1136565E-01
   0.8538573     -0.1027686E-01
   0.8707762     -0.9165672E-02
   0.8876519     -0.8034584E-02
   0.9044358     -0.6887851E-02
   0.9210328     -0.5733284E-02
   0.9372686     -0.4584932E-02
   0.9528497     -0.3466187E-02
   0.9673681     -0.2410030E-02
   0.9803942     -0.1452557E-02
   0.9916624     -0.6186220E-03
    1.000000      0.6798118E-10

Is there a solution for this problem?

John