with
  text_io,
  ada.integer_text_io,
  ada.long_float_text_io,
  ada.numerics,
  ada.numerics.long_elementary_functions;

procedure cylinder is

  type vec3 is array (1 .. 3) of long_float;

  ni,nj,nn,ne        : integer;
  radius,height,area : long_float;

  procedure scan (io : in text_io.file_type) is

    ascii : character;

  begin
    loop
      text_io.get(io,ascii);
      text_io.skip_line(io);
      exit when ascii = ':';
    end loop;
  end scan;

  function node (i,j : in integer) return integer is

  begin
    return ni * j + (i mod ni) + 1;
  end node;

  procedure vertex (i,j : in integer;
                    x   : out vec3) is

    use
      ada.numerics,
      ada.numerics.long_elementary_functions;

    t : long_float;

  begin
    t := 2.0 * pi / long_float(ni) * long_float(i mod ni);
    x(1) := radius * cos(t);
    x(2) := radius * sin(t);
    x(3) := (long_float(j) / long_float(nj) - 0.5) * height;
  end vertex;

  procedure initia is

    use
      text_io,
      ada.integer_text_io,
      ada.long_float_text_io;

    io : file_type;

  begin
    open(io, in_file, "cylinder.txt");

    -- The circumference is divided in ni parts.
    -- The height is divided in nj parts.
    -- The numbers ni and nj must be even.

    scan(io);
    get(io,radius);
    get(io,ni);
    get(io,height);
    get(io,nj);
    get(io,area);
    skip_line(io);

    close(io);

    -- nn, ne
    nn := ni * (nj + 1);
    ne := 3 * ni * nj + ni;
  end initia;

  procedure input is

    use
      text_io,
      ada.integer_text_io,
      ada.long_float_text_io;

    maxls,maxminor,maxlist,h    : integer;
    lstol,mntol,mnref,mu10,mu01 : long_float;
    e,glob,i,j                  : integer;
    x                           : vec3;
    io                          : file_type;

  begin
    create(io, out_file, "input.txt");

    -- solution parameters
    lstol := 1.0E-01;
    maxls := 20;
    mntol := 1.0E-05;
    mnref := 1.0;
    maxminor := 30 * nn;
    maxlist := 3;

    put(io,"        lstol     maxls             mntol             mnref  maxminor   maxlist");
    new_line(io);
    put(io,":");
    new_line(io);
    put(io,lstol,2,6,3);
    put(io,"     ");
    put(io,maxls,5);
    put(io,"     ");
    put(io,mntol,2,6,3);
    put(io,"     ");
    put(io,mnref,2,6,3);
    put(io,"     ");
    put(io,maxminor,5);
    put(io,"     ");
    put(io,maxlist,5);
    new_line(io);

    -- material parameters
    h := 2;
    mu10 := 0.3750;
    mu01 := -0.1250;

    new_line(io);
    new_line(io);
    put(io,"    h              mu10              mu01");
    new_line(io);
    put(io,":");
    new_line(io);
    put(io,h,5);
    put(io,"     ");
    put(io,mu10,2,6,3);
    put(io,"     ");
    put(io,mu01,2,6,3);
    new_line(io);

    -- node coordinates
    new_line(io);
    new_line(io);
    put(io," node           coord x           coord y           coord z");
    new_line(io);
    put(io,":");
    new_line(io);
    put(io,nn,5);
    new_line(io);

    new_line(io);

    for j in 0 .. nj loop
      for i in 0 .. ni - 1 loop
        vertex(i,j,x);
        put(io,node(i,j),5);
        put(io,"     ");
        put(io,x(1),2,6,3);
        put(io,"     ");
        put(io,x(2),2,6,3);
        put(io,"     ");
        put(io,x(3),2,6,3);
        new_line(io);
      end loop;
    end loop;

    -- elements
    new_line(io);
    new_line(io);
    put(io," elem      node      node              area");
    new_line(io);
    put(io,":");
    new_line(io);
    put(io,ne,5);
    new_line(io);

    e := 0;
    new_line(io);

    -- horizontal elements
    for j in 0 .. nj loop
      for i in 0 .. ni - 1 loop
        e := e + 1;
        put(io,e,5);
        put(io,"     ");
        put(io,node(i,j),5);
        put(io,"     ");
        put(io,node(i + 1,j),5);
        put(io,"     ");
        put(io,area,2,6,3);
        new_line(io);
      end loop;
    end loop;

    -- vertical elements
    for i in 0 .. ni - 1 loop
      for j in 0 .. nj - 1 loop
        e := e + 1;
        put(io,e,5);
        put(io,"     ");
        put(io,node(i,j),5);
        put(io,"     ");
        put(io,node(i,j + 1),5);
        put(io,"     ");
        put(io,area,2,6,3);
        new_line(io);
      end loop;
    end loop;

	-- diagonal elements
    j := 1;
    while j < nj loop

      i := 1;
      while i < ni loop
        e := e + 1;
        put(io,e,5);
        put(io,"     ");
        put(io,node(i,j),5);
        put(io,"     ");
        put(io,node(i + 1,j + 1),5);
        put(io,"     ");
        put(io,area,2,6,3);
        new_line(io);

        e := e + 1;
        put(io,e,5);
        put(io,"     ");
        put(io,node(i,j),5);
        put(io,"     ");
        put(io,node(i - 1,j + 1),5);
        put(io,"     ");
        put(io,area,2,6,3);
        new_line(io);

        e := e + 1;
        put(io,e,5);
        put(io,"     ");
        put(io,node(i,j),5);
        put(io,"     ");
        put(io,node(i - 1,j - 1),5);
        put(io,"     ");
        put(io,area,2,6,3);
        new_line(io);

        e := e + 1;
        put(io,e,5);
        put(io,"     ");
        put(io,node(i,j),5);
        put(io,"     ");
        put(io,node(i + 1,j - 1),5);
        put(io,"     ");
        put(io,area,2,6,3);
        new_line(io);

        i := i + 2;
      end loop;

      j := j + 2;
    end loop;

    -- known displacements
    new_line(io);
    new_line(io);
    put(io," node      axis             displ");
    new_line(io);
    put(io,":");
    new_line(io);
    put(io,6 * ni,5);
    new_line(io);

    new_line(io);

    for i in 0 .. ni - 1 loop
      glob := node(i,0);
      for axis in 1 .. 3 loop
        put(io,glob,5);
        put(io,"     ");
        put(io,axis,5);
        put(io,"      0.000000E+00");
        new_line(io);
      end loop;
    end loop;

    for i in 0 .. ni - 1 loop
      glob := node(i,nj);
      for axis in 1 .. 3 loop
        put(io,glob,5);
        put(io,"     ");
        put(io,axis,5);
        put(io,"      0.000000E+00");
        new_line(io);
      end loop;
    end loop;

    -- force
    new_line(io);
    new_line(io);
    put(io," node      axis             force");
    new_line(io);
    put(io,":");
    new_line(io);
    put(io,3 * ni * (nj - 1) - ni - 4 * (nj - 1),5);
    new_line(io);

    new_line(io);
    for j in 1 .. nj - 1 loop
      for i in 0 .. ni - 1 loop
        glob := node(i,j);
        vertex(i,j,x);
        for axis in 1 .. 3 loop
          if abs x(axis) > mntol then
            put(io,glob,5);
            put(io,"     ");
            put(io,axis,5);
            put(io,"     ");
            put(io,x(axis),2,6,3);
            new_line(io);
          end if;
        end loop;
      end loop;
    end loop;

    close(io);
  end input;

  procedure autocad is

    use
      text_io,
      ada.long_float_text_io;

    i,j : integer;
    x   : vec3;
    io  : file_type;

  begin
    create(io, out_file, "cylinder.scr");

    -- set running object snap mode to NONE
    put(io,"osmode 0");
    new_line(io);

    -- horizontal elements
    for j in 0 .. nj loop
      for i in 0 .. ni - 1 loop
        put(io,"3dpoly ");
        vertex(i,j,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        vertex(i + 1,j,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        new_line(io);
      end loop;
    end loop;

    -- vertical elements
    for i in 0 .. ni - 1 loop
      for j in 0 .. nj - 1 loop
        put(io,"3dpoly ");
        vertex(i,j,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        vertex(i,j + 1,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        new_line(io);
      end loop;
    end loop;

	-- diagonal elements
    j := 1;
    while j < nj loop

      i := 1;
      while i < ni loop
        put(io,"3dpoly ");
        vertex(i,j,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        vertex(i + 1,j + 1,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        new_line(io);

        put(io,"3dpoly ");
        vertex(i,j,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        vertex(i - 1,j + 1,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        new_line(io);

        put(io,"3dpoly ");
        vertex(i,j,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        vertex(i - 1,j - 1,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        new_line(io);

        put(io,"3dpoly ");
        vertex(i,j,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        vertex(i + 1,j - 1,x);
        put(io,x(1),1,8,3);
        put(io,",");
        put(io,x(2),1,8,3);
        put(io,",");
        put(io,x(3),1,8,3);
        put(io," ");
        new_line(io);

        i := i + 2;
      end loop;

      j := j + 2;
    end loop;

    put(io,"zoom ext ext");
    new_line(io);

    close(io);
  end autocad;

begin
  initia;
  input;
  autocad;
end cylinder;
