Вычислительный алгоритм построения численного конформного отображения четверти кольца на прямоугольник

#include <iostream>
#include <cmath>
 
using namespace std;
 
int main()
{
    int r1=1, r2=5;// variables for step 1
    int m=3, n=3;
    float y0j[n], xnj[n], xi0[m], yi0[m], xim[m], yim[m]; // variables for step 2
    float xij[m][n], yij[m][n]; //variables for step 3
    float xk[m][n], yk[m][n]; // variables for step 5
    double s[m][n],s_max; //variable for step 7
 
    //step 1
    float h1=3.14/(2*n);
    float h2=(log(r2/r1))/m;
 
    //step2 G3
    int x0j[m];
    int ynj[m];
 
    for(int j=0;j<n;j++)
     y0j[j]=r1+j*(r2-r1)/m;
    //step2 G4
    for(int j=0;j<n;j++)
     xnj[j]=r1+j*(r2-r1)/m;
    //step2 G1
    for(int i=0;i<m;i++)
    {
     xi0[i]=r1*sin(i*h1);
     yi0[i]=r1*cos(i*h1);
    }
     //step2 G2
    for(int i=0;i<m;i++)
    {
     xim[i]=r2*sin(i*h1);
     yim[i]=r2*cos(i*h1);
    }
    //step3
    for(int i=0;i<m;i++)
     for(int j=0;j<n;j++)
     {
      xij[i][j]=0.5*((xi0[i]+j*(xim[i]-xi0[i])/m)+(x0j[j]+i*(xnj[j]-x0j[j])/n));
      yij[i][j]=0.5*((yi0[i]+j*(yim[i]-yi0[i])/m)+(y0j[j]+i*(ynj[j]-y0j[j])/n));
     }
     //step 5
     float psi=h1/h2;
     for(int i=1;i<m-1;i++)
      for(int j=1;j<n-1;j++)
      {
       xk[i][j]=(xij[i-1][j]+xij[i+1][j]+pow(psi,2)*((xij[i][j-1]+xij[i][j+1])))/(2*1+pow(psi,2));
       yk[i][j]=(yij[i-1][j]+yij[i+1][j]+pow(psi,2)*((yij[i][j-1]+yij[i][j+1])))/(2*1+pow(psi,2));
      }
      //step 6
      float psi_temp;
      for(int i=0;i<m-1;i++)
       for(int j=0;j<n-1;j++)
       {
        psi_temp=(sqrt(pow((xij[i+1][j]-xij[i][j]),2)+pow((yij[i+1][j]-yij[i][j]),2))+sqrt(pow((xij[i+1][j+1]-xij[i][j+1]),2)+pow((yij[i+1][j+1]-yij[i][j+1]),2)))/(sqrt(pow((xij[i+1][j]-xij[i][j]),2)+pow((yij[i+1][j]-yij[i][j]),2))+sqrt(pow((xij[i+1][j+1]-xij[i][j+1]),2)+pow((yij[i+1][j+1]-yij[i][j+1]),2)));
       }
       psi=psi_temp*(1/n*m);
      //step 7
      for(int i=1;i<m;i++)
       for(int j=1;j<n;j++)
       {
        s[i-1][j-1]=sqrt(pow((xk[i][j]-xk[i-1][j-1]),2)+(pow((yk[i][j]-yk[i-1][j-1]),2)));
       }
       s_max=s[0][0];
       for(int i=0;i<m;i++)
        for(int j=0;j<n;j++)
        if(s[i][j]>s_max)
         s_max=s[i][j];
      //results
      cout<<"Step 1:"<<endl;
      cout<<"h1="<<h1<<endl;
      cout<<"h2="<<h2<<endl;
      cout<<"Step 2:"<<endl;
      for(int j=0;j<n;j++)
       cout<<y0j[j]<<endl;
      for(int j=0;j<n;j++)
       cout<<xnj[j]<<endl;
      for(int i=0;i<m;i++)
      {
       cout<<xi0[i]<<endl;
       cout<<yi0[i]<<endl;
      }
      for(int i=0;i<m;i++)
      {
       cout<<xim[i]<<endl;
       cout<<yim[i]<<endl;
      }
      cout<<"Step 3:"<<endl;
      for(int i=0;i<m;i++)
       for(int j=0;j<n;j++)
       {
        cout<<xij[i][j]<<endl;
        cout<<yij[i][j]<<endl;
       }
       cout<<"Step 5:"<<endl;
       cout<<psi<<endl;
       for(int i=0;i<m;i++)
        for(int j=0;j<n;j++)
        {
         cout<<xk[i][j]<<endl;
         cout<<yk[i][j]<<endl;
        }
        cout<<"s="<<s_max;
    return 0;
}