%
% Program proj_symprjs.m extracts the projection matrices from SYMPRJS binary
% file output.
%
% Version 1.0, 2006-08-06
%
in_proj = fopen('wurtzite_proj','r');
out_proj = fopen('symprjs_proj_wurtzite','a');
no_wave_vectors = fread(in_proj,1,'int8');
nel = fread(in_proj,1,'int8');
for I = 1:nel
  nat(I) = fread(in_proj,1,'int8');
end
for I = 1:nel
  lmax(I) = fread(in_proj,1,'int8');
end
matrix_order = 0;
for I = 1:nel
  g = 0;
  for J = 1:(lmax(I) + 1)
    g = g + (2*(J - 1) + 1);
  end
  matrix_order = nat(I)*g + matrix_order;
end
matrix_order = 2*matrix_order;
fwrite(out_proj, no_wave_vectors, 'int8');
fwrite(out_proj, matrix_order, 'int16');
column_index(1,1) = 1;
row_index(1,1) = 1;
col = 1;
for I = 1:nel
  J = nat(I);
  for K = 2:(lmax(I) + 1)
    l = 2*J*(2*(K - 2) + 1);
    col = col + l;
    column_index(K,I) = col;
    row_index(K,I) = col;
  end
  if I < nel
    col = col + 2*J*(2*lmax(I) + 1);
    column_index(1,I+1) = col;
    row_index(1,I+1) = col;
  end
end 
wave_vector = 0;
while wave_vector < no_wave_vectors
  wave_vector = wave_vector + 1;
  tmat = zeros(matrix_order);
  for I = 1:nel
    for II = 1:(lmax(I)+1)
      column_index_temp(II,I) = column_index(II,I);
    end
  end
  while column_index_temp(lmax(nel)+1,nel) < matrix_order 
    ichem = fread(in_proj,1,'int8');
    L = fread(in_proj,1,'int8');
    itotal = fread(in_proj,1,'int8');
    ndi = fread(in_proj,1,'int8');
    for I = 1:itotal
      for II = 1:ndi
        for III = 1:2
          tmatridel(III) = fread(in_proj,1,'float64');
        end
        tmat(II+row_index(L+1,ichem)-1,column_index_temp(L+1,ichem)) = tmatridel(1) + i*tmatridel(2);
      end
      column_index_temp(L+1,ichem) = column_index_temp(L+1,ichem) + 1;
    end
  end
  tmatriorg(1,1:matrix_order,1:matrix_order) = real(tmat(1:matrix_order,1:matrix_order));
  tmatriorg(2,1:matrix_order,1:matrix_order) = imag(tmat(1:matrix_order,1:matrix_order));
  fwrite(out_proj, tmatriorg(1:2,1:matrix_order,1:matrix_order), 'float64');
end
fclose(in_proj);
fclose(out_proj);
fprintf('proj_symprjs has finished executing');