%
% Program prodats prepares datafiles to be used by program SYMPRJS/SYMPWS.
%
% Version 1.0, 060423
%
fid = fopen('prodats.dat','r');
fid2 = fopen('prodatsres.m','a');
Oh = fscanf(fid,'%g %g %g',[3,24]);
for I = 1:24
  Oh(4,I) = 0;
  Oh(1:3,I) = Oh(1:3,I)*pi;
  Oh(1:3,I+48) = Oh(1:3,I);
  Oh(4,I+48) = 1;
end
Oh(1:4,25:48) = Oh(1:4,1:24);
Oh(1:4,73:96) = Oh(1:4,49:72);
fprintf(fid2, 'Rotation angles for double group Oh\n\n');  
for I = 1:96
  fprintf(fid2, '%2d    %7.4f %7.4f %7.4f %7.4f\n', I, Oh(1:4,I));
end
fprintf(fid2, '\n');
D6h = fscanf(fid,'%g %g %g',[3,12]);
for I = 1:12
  D6h(4,I) = 0;
  D6h(1:3,I) = (D6h(1:3,I)/3)*pi;
  D6h(1:3,I+24) = D6h(1:3,I);
  D6h(4,I+24) = 1;
end
D6h(1:4,13:24) = D6h(1:4,1:12);
D6h(1:4,37:48) = D6h(1:4,25:36);
fprintf(fid2, '\n');
fprintf(fid2, 'Rotation angles for double group D6h\n\n');
for I = 1:48
  fprintf(fid2, '%2d    %7.4f %7.4f %7.4f %7.4f\n', I, D6h(1:4,I));
end
fprintf(fid2, '\n');
MOh = fscanf(fid,'%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d',[48,48]);
MOh = MOh.';
MOh(1:48,49:96) = MOh(1:48,1:48) + 48;
MOh(49:96,1:48) = MOh(1:48,1:48) + 48;
MOh(49:96,49:96) = MOh(1:48,1:48);
MD6h = fscanf(fid,'%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d',[24,24]);
MD6h = MD6h.';
MD6h(1:24,25:48) = MD6h(1:24,1:24) + 24;
MD6h(25:48,1:24) = MD6h(1:24,1:24) + 24;
MD6h(25:48,25:48) = MD6h(1:24,1:24);
fprintf(fid2, '\n');  
fprintf(fid2, 'The double group Oh multiplication table:\n\n');
for I = 1:96
  fprintf(fid2, '%2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d\n', MOh(I,1:96));
  fprintf(fid2, '\n');
end
fprintf(fid2, '\n\n');
fprintf(fid2, 'The double group D6h multiplication table:\n\n');
for I = 1:48
  fprintf(fid2, '%2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d\n', MD6h(I,1:48));
  fprintf(fid2, '\n');
end
fprintf(fid2, '\n');
npgo = fscanf(fid,'%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d', [36,2]);
npgo = npgo.';
for I = 1:2
  fprintf(fid2,'npgo(%1d,1:36)\n\n', I);
  fprintf(fid2,'%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n',npgo(I,1:36));
  fprintf(fid2, '\n\n');
end
fprintf(fid2,'The group elements of the 36 crystallographic point groups:\n\n');
for I = 1:36
  K = npgo(1,I);
  L = npgo(2,I);
  L2 = L + K - 1;
  nge(L:L2) = fscanf(fid,'%d %d %d %d %d %d %d %d', K);
  fprintf(fid2, '%2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d\n', nge(L:L2));
  fprintf(fid2, '\n');
end
fprintf(fid2, '\n');
for I = 1:24
  fi = Oh(1,I);
  theta = Oh(2,I);
  psi = Oh(3,I);
  fprintf(fid2, '\n');
  fprintf(fid2,'%3d ANGLES %7.5f %7.5f %7.5f\n', I, fi, theta, psi);
  for K = 1:4
    L = K - 1;
    fclose(fid2);
    D = dmatr(L,fi,theta,psi);
    fid2 = fopen('prodatsres.m', 'a');
    N = 2*L + 1;
    K1 = 1;
    if L ~= 1
      if L == 2
        K1 = 11;
      end
      if L == 3
        K1 = 36;
      end
    else
      rcgr3(1:3,1:3,I) = D(1:3,1:3);
      rcgr3(1:3,1:3,I+24) = -D(1:3,1:3);
      K1 = 2;
    end
    for L1 = 1:N
      for L2 = 1:N
        K2 = K1 + (L1 - 1)*N + L2 - 1;
        ldrmm(I,K2) = D(L1,L2);
        ldrmm(I+48,K2) = D(L1,L2)*(-1)^L;
      end
    end
  end
end
ldrmm(25:48,1:84) = ldrmm(1:24,1:84);
ldrmm(73:96,1:84) = ldrmm(49:72,1:84);
for I = 1:12
  fi = D6h(1,I);
  theta = D6h(2,I);
  psi = D6h(3,I);
  fprintf(fid2, '\n');
  fprintf(fid2,'%3d ANGLES %7.5f %7.5f %7.5f\n', I, fi, theta, psi);
  for K = 1:4
    L = K - 1;
    fclose(fid2);
    D = dmatr(L,fi,theta,psi);
    fid2 = fopen('prodatsres.m','a');
    N = 2*L + 1;
    K1 = 1;
    if L ~= 1
      if L == 2
        K1 = 11;
      end
      if L == 3
        K1 = 36;
      end
    else
      rcgr3(1:3,1:3,I+48) = D(1:3,1:3);
      rcgr3(1:3,1:3,I+60) = -D(1:3,1:3);
      K1 = 2;
    end
    for L1 = 1:N
      for L2 = 1:N
        K2 = K1 + (L1 - 1)*N + L2 - 1;
        ldrmm(I+96,K2) = D(L1,L2);
        ldrmm(I+120,K2) = D(L1,L2)*(-1)^L;
      end
    end
  end
end
ldrmm(109:120,1:84) = ldrmm(97:108,1:84);
ldrmm(133:144,1:84) = ldrmm(121:132,1:84);
irt2 = 1/sqrt(2);
Q1 = [irt2 0 -irt2; -i*irt2 0 -i*irt2; 0 1 0];
Q = Q1';
knr2 = 24;
knr = 0;
for I = 1:72
  if I > 24
    knr = 24;
  end
  if I > 48
    knr = 48;
    knr2 = 12;
  end
  if I > 60
    knr = 60;
  end
  res(1:3,1:3) = Q1(1:3,1:3)*rcgr3(1:3,1:3,I);
  res2(1:3,1:3) = res(1:3,1:3)*Q(1:3,1:3);
  rgr3(1:3,1:3,I+knr) = res2(1:3,1:3);
  rgr3(1:3,1:3,I+knr+knr2) = res2(1:3,1:3);
end
fprintf(fid2,'\n\n');
fprintf(fid2,'Rotation/inversion matrices\n\n');
for I = 1:144
  fprintf(fid2,'%3d\n', I);
  for K1 =1:3
    fprintf(fid2,'%7.4f %7.4f %7.4f\n', rgr3(K1,1:3,I));
  end
  fprintf(fid2,'\n');
end
NP = primen;
fprintf(fid2, '\n\n');
fprintf(fid2, 'One hundred prime numbers are calculated:\n\n');
fprintf(fid2, '%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', NP);
fprintf(fid2, '\n');
for I = 1:24
  fi = Oh(1,I);
  theta = Oh(2,I);
  psi = Oh(3,I);
  spinm(2,2,I) = cos(theta/2)*exp(i*(fi+psi)/2);
  spinm(1,1,I) = conj(spinm(2,2,I));
  spinm(2,1,I) = sin(theta/2)*exp(i*(fi-psi)/2);
  spinm(1,2,I) = -conj(spinm(2,1,I));
end
spinm(1:2,1:2,25:48) = -spinm(1:2,1:2,1:24);
spinm(1:2,1:2,49:72) = spinm(1:2,1:2,1:24);
spinm(1:2,1:2,73:96) = -spinm(1:2,1:2,1:24);
for I = 1:12
  fi = D6h(1,I);
  theta = D6h(2,I);
  psi = D6h(3,I);
  spinm(2,2,I+96) = cos(theta/2)*exp(i*(fi+psi)/2);
  spinm(1,1,I+96) = conj(spinm(2,2,I+96));
  spinm(2,1,I+96) = sin(theta/2)*exp(i*(fi-psi)/2);
  spinm(1,2,I+96) = -conj(spinm(2,1,I+96));
end
spinm(1:2,1:2,109:120) = -spinm(1:2,1:2,97:108);
spinm(1:2,1:2,121:132) = spinm(1:2,1:2,97:108);
spinm(1:2,1:2,133:144) = -spinm(1:2,1:2,97:108);
fwritergr3 = fopen('rgr3','a');
fwriteldrmm = fopen('ldrmm','a');
fwritespinm = fopen('spinm','a');
fwritenpgo = fopen('npgo','a');
fwritenge = fopen('nge','a');
fwriteMOh = fopen('MOh','a');
fwriteMD6h = fopen('MD6h','a');
fwriteNP = fopen('NP','a');
for I = 1:144
  for K1 =1:3
    for K2 = 1:3  
      fwrite(fwritergr3,rgr3(K1,K2,I), 'float64');
    end  
  end
end
ldrmmorg(1,1:144,1:84) = real(ldrmm(1:144,1:84));
ldrmmorg(2,1:144,1:84) = imag(ldrmm(1:144,1:84));
for I = 1:144
  for K1 = 1:84
    for K2 = 1:2
      fwrite(fwriteldrmm, ldrmmorg(K2,I,K1),'float64');
    end
  end
end
spinmorg(1,1:2,1:2,1:144) = real(spinm(1:2,1:2,1:144));
spinmorg(2,1:2,1:2,1:144) = imag(spinm(1:2,1:2,1:144));
for I = 1:144
  for K1 = 1:2
    for K2 =1:2  
      for K3 = 1:2
        fwrite(fwritespinm,spinmorg(K3,K1,K2,I), 'float64');
      end  
    end
  end
end  
for I = 1:2
  fwrite(fwritenpgo,npgo(I,1:36), 'int16');
end
fwrite(fwritenge,nge(1:736), 'int16');
for I = 1:96
  fwrite(fwriteMOh,MOh(I,1:96), 'int8'); 
end
for I = 1:48
  fwrite(fwriteMD6h,MD6h(I,1:48), 'int8'); 
end
fwrite(fwriteNP,NP(1:100), 'int16');
fclose(fwritergr3);
fclose(fwriteldrmm);
fclose(fwritespinm);
fclose(fwritenpgo);
fclose(fwritenge);
fclose(fwriteMOh);
fclose(fwriteMD6h);
fclose(fwriteNP);
fclose(fid);
fclose(fid2);
fprintf('Prodats has finished executing');