Вычисление площади кривой

clear
close
x=[77 80 84 85 86 90 91.5 93.5 92 90 87 80 75 70 65 60 55 50 45 35.5 30 25 23 20 19 15 10 5 1 -5 -10 -15 -20 -25 -30 -34 -40 -44 -50 -55 -60 -65 -70 -75 -80 -85 -90 -94 -97.5 -99 -100.5 -101 -101.5 -101 -100 -98 -96 -90 -86.5 -83 -79 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -19.5 -10 -5 4 10 15 20 30 40 50 55 60.5 65 70 77];
y=[0 5 10 15 20 30 35 40 50 55 60 67.5 71 74 75 75.5 76.5 76 74 70 66 60 55 50 46 43.5 43 42 42 41.8 42 42.5 45 50 55 60 65 70 76 77 78.8 79.5 80 78 76 74 70 65 60 56 50 45 40 35 35 30 25 18 15 10 6 -7 -11 -15.5 -18 -21 -23 -25 -26 -27.5 -28 -29 -29.5 -30 -30 -30 -29.8 -29.5 -27.5 -25 -20 -17 -14 -10 -7 0];
phi = atan (y./x);
for I=1: length(x)
if x(I)<0
    phi(I) = phi(I) + pi;
    else
        if y(I)<0
            phi(I) = phi(I)  + 2*pi;
        end
    end
end
             phi(length (x)) = 2*pi;
           
          
               phin=0:pi./180:2*pi;
        xnl=interp1(phi,x,phin,'linear');
         ynl=interp1(phi,y,phin,'linear');
         xnc=interp1(phi,x,phin,'cubic');
         ync=interp1(phi,y,phin,'cubic');
         xns=interp1(phi,x,phin,'spline');
         yns=interp1(phi,y,phin,'spline');
         
         subplot(3,2,1);
            axis('off');text(0,1,'Grafik, Gr.5131');
            subplot(3,2,[2 6]);
         plot(x,y,xnl,ynl,xnc,ync,xns,yns);
         axis('equal');
         grid;
         subplot(3,2,3);
         plot(phi,x,phin,xnl,phin,xnc,phin,xnl);
               grid;
         subplot(3,2,5);
         plot(phi,y,phin,ynl,phin,ync,phin,yns);
         grid;