⟦37004999f⟧

TextFile

;gosav
lookup struktur
if ok.no
(r=algol index.no
rename r.struktur
permanent struktur.17)
\f





STRUKTURPROGRAM 17-5-1974. GOS.

begin
comment The input rules are:
      1) A textstring in < >.
      2) The number of atoms in the structure (N).
      3) For each atom: Up to three coordinates (last del.: ,)
          and the appropiate uncertainties (last del.: nl.).
               x1  y1  z1,        dx1  dy1  dz1
               x2  .  .  .
      4) Any number of the following types of constructions:
               p  dist:  i1 j1   i2 j2   .  .  ip jp
               p  angle: i1 j1 k1   .  .  .  ip jp kp
         or    p  outo:  i1 j1 k1 l1 .  .  ip jp kp lp
         indicating the desired bondlengths and angles.
      5) The input is terminated by p = 0;
integer i,j,k,l,m,p,q,N;
boolean nl, sp, closeres;
array H(1:1), head(1:12);
zone res(128,1,stderror);
closeres:= outmedium(res);
nl:= false add 10; sp:= false add 32;

new_data:
readhead(in,head,1); read(in,N);
i:= 1;   write(res,<:<12>:>,nl,4,string head(increase(i)),nl,2,<:
CARTESIAN COORDINATES (Å)::>,nl,2);

begin
real type, rad;
array X, dX(1:N,1:3);
real procedure vector(e,i,j);
   value i,j; integer i,j; array e;
   begin
   real l,r; integer k;
   r:= 0;
   for k:=1 step 1 until 3 do begin
      e(k):= l:= X(j,k)-X(i,k); r:= r+ l*l end;
   r:= vector:= sqrt(r); r:= 1/r;
   for k:= 1 step 1 until 3 do e(k):= e(k)*r
end vector;
real procedure dot(a,b);
   array a,b;
   begin real d; integer k;
   d:= 0; for k:=1 step 1 until 3 do d:= d+ a(k)*b(k);
   dot:= d
end dot;
real procedure cross(a,b,c);
   array a,b,c;
   begin real r; integer k;
   c(1):= a(2)*b(3) - a(3)*b(2);
   c(2):= a(3)*b(1) - a(1)*b(3);
   c(3):= a(1)*b(2) - a(2)*b(1);
   r:= 0; for k:=1 step 1 until 3 do r:= r+ c(k)**2;
   cross:= r
end cross;
rad:= 180/pi;  k:= 0;
for i:=1 step 1 until N do begin
   for j:=1,j+1 while k<>44 do begin
      read(in,X(i,j)); repeatchar(in); readchar(in,k) end;
   for j:=j step 1 until 3 do X(i,j):= 0;
   for j:=1,j+1 while k<>10 do begin
      read(in,dX(i,j)); repeatchar(in); readchar(in,k) end;
   for k:=1 step 1 until 3 do if X(i,k)=0
   then write(res,<:      0       :>)
   else write(res,<<____-dd.dddddd>,X(i,k));
   write(res,nl,1);
   for k:=1 step 1 until j-1 do
        write(res,<<______d.d00000>,dX(i,k));
   write(res,nl,2);
   for j:=j step 1 until 3 do dX(i,j):= 0;
   for j:= 1 step 1 until 3 do dX(i,j):= dX(i,j)**2
end i;
read(in,p);
newtype:
readstring(in,H,1); type:= H(1);

if type = real <:dist:> then begin
real r, dr;
array e(1:3);
for m:=1 step 1 until p do begin
   read(in,i,j);  r:= vector(e,i,j);  dr:= 0;
   for q:=1 step 1 until 3 do dr:= dr+ e(q)**2*(dX(i,q)+dX(j,q));
   write(res,nl,1,<:
distance :>,<<dd>,i);
   if j>9 then write(res,<<-dd>,-j) else
   write(res,<<-d>,-j,sp,1);
   write(res,<: = :>,<<d.ddddd>,r,sp,2,<<d.d0000>,sqrt(dr))
end;  type:= 0;  write(res,nl,1)
end;

if type = real <:angle:> then begin
real cos,sin,ri,rk,v,dv,a,b;
array ei,ek(1:3);
for m:=1 step 1 until p do begin
   read(in,i,j,k);
   ri:= vector(ei,j,i); rk:= vector(ek,j,k); cos:= dot(ei,ek);
   sin:= sqrt(1-cos*cos);  v:= arctan(sin/cos);
   if v<0 then v:= v+pi; ri:= 1/(sin*ri); rk:= 1/(sin*rk);  dv:= 0;
   for q:=1 step 1 until 3 do begin
      a:= (ei(q)*cos- ek(q))*ri;
      b:= (ek(q)*cos- ei(q))*rk;  sin:= (a+b)**2;
      dv:= dv+ a*a*dX(i,q) + b*b*dX(k,q) + sin*dX(j,q)  end;
   write(res,nl,1,<:
angle :>,<<dd>,i);  q:= 0;
   if j>9 then write(res,<<-dd>,-j) else
   begin write(res,<<-d>,-j); q:= q+1 end;
   if k>9 then write(res,<<-dd>,-k) else
   begin write(res,<<-d>,-k); q:= q+1 end;
   write(res,sp,q,<: = :>,<<ddd.dd>,v*rad,sp,2,
   <<d.d0>,sqrt(dv)*rad,<: degrees.:>)
end;  type:= 0;  write(res,nl,1)
end;

if type = real <:outo:> then begin
real ri,rj,rl,sin,csc,cot,a,b,c,d,v,dv;
array ei,ej,el,cij,cjl,cli(1:3);
for m:=1 step 1 until p do begin
   read(in,i,j,k,l);
   ri:= 1/vector(ei,k,i);
   rj:= 1/vector(ej,k,j);
   rl:= 1/vector(el,k,l);
   csc:= 1/sqrt(cross(ei,ej,cij));
   cross(ej,el,cjl);  cross(el,ei,cli);
   cot:= dot(ei,ej)*csc;  sin:= dot(cij,el)*csc;
   a:= 1/sqrt(1-sin*sin);  v:= arctan(a*sin);
   dv:= 0; ri:= ri*a*csc;  rj:= rj*a*csc;  rl:= rl*a;
   for q:=1 step 1 until 3 do begin
      a:= (cjl(q)-(ei(q)*csc-ej(q)*cot)*sin)*ri;
      b:= (cli(q)-(ej(q)*csc-ei(q)*cot)*sin)*rj;
      c:= (cij(q)*csc-el(q)*sin)*rl;
      d:= -a-b-c;
      dv:= dv + a*a*dX(i,q)+b*b*dX(j,q)+c*c*dX(l,q)+d*d*dX(k,q);
   end;
   write(res,nl,1,<:
outo :>,<<dd>,i); q:= 0;
   if j>9 then write(res,<<-dd>,-j) else
   begin write(res,<<-d>,-j); q:= q+1 end;
   if k>9 then write(res,<<-dd>,-k) else
   begin write(res,<<-d>,-k); q:= q+1 end;
   if l>9 then write(res,<<-dd>,-l) else
   begin write(res,<<-d>,-l); q:= q+1 end;
   write(res,sp,q,<: = :>,<<-ddd.dd>,v*rad,sp,2,
   <<d.d0>,sqrt(dv)*rad,<: degrees.:>)
end;  type:= 0;  write(res,nl,1)
end;

comment additional types can be included here;

if type<>0 then begin
write(res,nl,2,string type,<: : Type undefined.:>);
goto stop
end;
read(in,p); if p<>0 then goto newtype
end X;
stop:
write(res,nl,5,<:<25>:>); close(res,closeres)
end
▶EOF◀
DataMuseum.dk

RC4000/8000/9000

⟦37004999f⟧ TextFile

Derivation

TextFile
