The following ALGOL program of J.P. Dahl is presented here 
as an example of how programming was done around 1965.


Please note: for simplicity of description 
the original  was replaced by →, 
and the original  consequently by ←.





begin comment
        Atomic Program no. 1        September 28, 1965 (J. P. Dahl)
                                    June 20, 1967
            Generates radial charge density from a set of analytic atomic
            orbitals on a uniform or a Herman-Skillman mesh;
         
integer Zatom, number of basis sets, basis set number, M, N, norm,
        P2place, limit, dk, mesh type;
real mu, r, delR;
         
         
real procedure Rvalue (n);
        integer n;
        begin integer block,n1;
         if mesh type=1 then 
           begin 
              block := entier((n-1)/40);
              nl :=2↑block;
              Rvalue := mu*(0.1*(n1-1) + (n-block*40-1)*n1*0.0025)
            end
            else
            Rvalue :=(n-1)*delR
         end;
         
   outtext(→< Radial charge density [ ←);
   outcopy(→<[ ]←);
   outtext(→<]←);
   input(Zatom, norm, mesh type, delR);
   writetext(→< How many basis sets: ←);
   number of basis sets := typein;
   writetext(→< Specify outputmesh (dk=1,2,3 or 4): ←);
   dk := typein;
         
   mu := 0.88534138/Zatom↑(1/3);
         
begin integer k; array P2[1:441];
       for k:= 1 step dk until 441 do P2[k] := 0;
          P2place := drumplace;
          to drum (p2)
end;


   basis set number := 0;
         
  A:

    basis set number := basis set number + 1;
    if basis set number > number of basis sets then go to Punch;
         M := inone
         write(→ndd←), writetext(→< Basis set number←), basis set number);
         writetext(→< How many orbitals:←); 
         N := typein;
         
         limit := (M*(M+1))/2;
            begin integer s,t,u;
            integer array n[1:M];
            real array z, NN[1:M], S[1:limit];



            integer procedure Factorial(m); integer m;
               begin integer F;
                       F := if m=1 then 1 else m*Factorial(m-1);
                       Factorial := F
             end;
         

          for s := 1 step 1 until M do 
                 begin
                  input(n[s],z(s]);
                  NN[s] := sqrt((2*z(s])↑(2*n[s]+1) / Factorial(2*n[s]))
                  end;

          u := 0;

          for s := 1 step 1 until M do
          for t := s step 1 until M do 
             begin 
                 u := u+1;
                 S[u] := Factorial(n[s]+n[t])/(z[s] + z[t])↑(n[s]+n[t]+1);
                 if norm=0 then S[u]:= NN[s]*NN[t]*S[u]
              end;



            begin integer i, factor; real dsum, sum;
            array C[1:N,1:M], nn[1:N];
            input(C);
 
            writetext(→< Please specify number of electrons in each orbital 
                                  (real type numbers): ←);
           for i:= 1 step 1 until N do
              nn[i] := typein;
              writetext(→< Thank you. ←);
         
          for i := 1 step 1 until N do
                  begin
                  u := 0;
                  sum := 0;
                  for s := 1 step 1 until M do 
                 for t := s step 1 until M do 
                 begin 
                  u := u+1;
                  factor := if s=t then 1 else 2;
                  sum:= sum + C[i,s] * C[i,t] * S[u] * factor
                  end;

                  write(→dd←), writetext(→<
                  Orbital number←), i, writetext(→< Normalization : ←),
                  write(→nd.dddd dddd←), sum))
                  end;

            u := 0;

            for s := 1 step 1 until M do
            for t := s step 1 until M do
                 begin 
                 u := u+1;
                 sum := 0;
                 factor := if s=t then 1 else 2;
                 for i:=1 step 1 until N do
                      begin 
                      dsum := nn[i]*C[i,t]*C[i,s]*factor;
                      if norm=0 then dsum := dsum*NN[s]*NN[t];
                      sum := sum + dsum
                  end;
            S[u] := sum
           end
       end;
         

  begin integer k; real sum; array P2[1: 441];
       drumplace := P2place;
       from drum (P2);
       for k:=1 step dk until 441 do 
                 begin 
                 r := Rvalue(k);
                 u :=0; sum := 0;
                  for s :=1 step 1 until M do
                  for t :=s step 1 until M do 
                     begin   
                     u := u+1;
                     sum := sum + S[u]* r↑(n[s] + n[t]) * exp((-z[s] - z[t] ) * r)
                  end;
                  P2[k] := P2[k] + sum
            end;
         
            drumplace := P2place;
            to drum (P2)
            end 
   end;
         
            go to A;
         

   Punch:

   begin integer k; array [1:441];
         outclear;
         outtext(→< Meshpoint Radius (a u) Radial charge density ←);
         drumplace := P2place;
         from drum (P2);
         for k:= 1 step dk until 441 do
            begin
            r := Rvalue(k);
            outcr;

            output(→nddddd←, k); outtext(→<,←); outsp(4);
            output(→ndd.dddd dddd←, r, outtext(→<,←),  outsp(5), P2[k])
            end;
         outsum
         end
         end;