如何解决以下报错

当我运行程序时出现下列报错，然而我在power函数中使用的nihejieguo为拟合产生，符合power函数的调用规则，现在如何解决？

Index exceeds the number of array elements. Index must not exceed 1.

出错 solve11>power (第 60 行)

A=x(1).*(x(2).*xdata.^4+x(3).*xdata.^3+x(4).*x.^2+x(5).*x+x(6))

clear

clc

%使用lsqcurvefit函数进行太阳光在大气层外表面的辐射强度拟合

data1=readtable("Afujian1.xlsx");

data2=readtable("Afujian2.xlsx");

phiwei=2*pi*(30.583/360);

phijin=2*pi*(114.317/360); %导入数据

x01=[1353,0]; %设定迭代初值

xdata1=data2.month*30-15;

ydata1=data2.I0_W_m2_;

x1=lsqcurvefit(@sunpower,x01,xdata1,ydata1); %求解

figure( 'Name', '太阳光在大气层外表面辐射强度' );

plot([1:365],sunpower(x1,[1:365]));

hold on;

scatter(xdata1,ydata1);

title('太阳光在大气层外表面辐射强度');

xlabel('日期序号');

ylabel('I_0 W/m^2')

Iavg=sum(sunpower(x1,[1:365]))./365; %平均辐射强度

st=data1.x___h;

omiga=30.*(st-12).*pi./360;

N=143;

h=1000;

delta=0.39795.*cos(0.98563.*(N-173)); %太阳赤纬

xdata2=data1.x___h;

sinH=sin(phiwei).*sin(delta)+cos(phiwei).*cos(delta).*cos(omiga);

cosH=sqrt(1-sinH.*sinH);

const=sunpower(x1,143);

xdata2=st

ydata2=data1.x_______W_m2;

[fitresult, gof] = createFit(xdata2, ydata2);

nihejieguo=coeffvalues(fitresult);

N0=[15,46,74,105,135,166,196,227,258,288,319,349]';

t=[8:0.05:17];

I0=sunpower(x1,N0);

delta0=0.39795.*cos(0.98563.*(N0-173));

omiga0=30.*(t-12).*pi./360;

sinH0=sin(phiwei)*sin(delta0)+cos(phiwei)*cos(delta0)*cos(omiga0);

cosH0=sqrt(1-sinH0.*sinH0);

cosG0=(sin(delta0)-sinH0.*sin(phiwei))./(cosH0.*cos(phiwei));

sinG0=sqrt(1-cosG0.*cosG0);

theta=0;

phi1=pi./9;

yuxuansunshi1=abs(cosH0.*cosG0.*cos(phi1).*cos(theta)+cosH0.*sinG0.*cos(phi1).*sin(theta)+sin(phi1).*sinH0);

I1=I0.*yuxuansunshi1.*power(nihejieguo,t)./const;

phi2=2*pi./9;

yuxuansunshi2=abs(cosH0.*cosG0.*cos(phi2).*cos(theta)+cosH0.*sinG0.*cos(phi2).*sin(theta)+sin(phi2).*sinH0);

phi3=pi./3;

yuxuansunshi3=abs(cosH0.*cosG0.*cos(phi3).*cos(theta)+cosH0.*sinG0.*cos(phi3).*sin(theta)+sin(phi3).*sinH0);

function F=sunpower(x,xdata)

F=x(1)*(1-0.0167*cos(0.0172*xdata+x(2))).*(1-0.0167*cos(0.0172*xdata+x(2)))

end

%太阳发出热辐射到地球大气层外表面的单位面积辐射量

function A=power(x,xdata)

A=x(1).*(x(2).*xdata.^4+x(3).*xdata.^3+x(4).*x.^2+x(5).*x+x(6))

end

%下列函数为curve fitting toolbox 自动生成代码

function [fitresult, gof] = createFit(st, ydata2)

[xData, yData] = prepareCurveData( st, ydata2 );

% 设置 fittype 和选项。

ft = fittype( 'a.*(p1*x^4+p2*x^3+p3*x^2+p4*x+p5)', 'independent', 'x', 'dependent', 'y' );

opts = fitoptions( 'Method', 'NonlinearLeastSquares' );

opts.Display = 'Off';

opts.Lower = [1200 -Inf -Inf -Inf -Inf -Inf];

opts.StartPoint = [1300 0.91719366382981 0.757200229110721 0.753729094278495 0.380445846975357 0.567821640725221];

opts.Upper = [2600 Inf Inf Inf Inf Inf];

% 对数据进行模型拟合。

[fitresult, gof] = fit( xData, yData, ft, opts );

% 绘制数据拟合图。

figure( 'Name', '地表太阳辐射强度' );

h = plot( fitresult, xData, yData );

legend( h, 'ydata2 vs. st', '拟合结果', 'Location', 'NorthEast', 'Interpreter', 'none' );

% 为坐标区加标签

xlabel( 'st', 'Interpreter', 'none' );

ylabel( 'ydata2', 'Interpreter', 'none' );

title('地表太阳辐射强度');

grid on