function write_covar2pipe(C, C_axes,fname_out)
% write_covar2pipe(C, C_axes,fname_out)
% write results of covariance to an NMRPipe format file named fname_out
% covar_out is result (either struct or string naming a temp file) produce as output of covar.m

if(nargin ~= 3)
    help write_covar2pipe
    return
end

num_dims = sum(C_axes(26,:) > 1);
load template_pipeheader; % header read from a typical NMR Pipe file and stored for use here

header_out = reset_header_data(C_axes, template_pipeheader, num_dims);

% open output file and write NMR pipe header

pipe_outfile = fopen(fname_out,'w');
fwrite(pipe_outfile,header_out,'float32');

if(isstr(C)) % output from covar rather than being covar spec is string pointing to file with output info
    [S_size U_size V_size C_size svd_power Z_flag sigma_A sigma_B temp_array rows2use cols2use S U V] = read_tmp_file(C);
    
    if(Z_flag ~= 2) % we only need V in case of mrs
        clear V;
    end

	clear V_size; % we won't actually use these
    clear temp_array;
    
    S = S.^(svd_power/2); % raise S to appropriate power for SVD

	rows2use(1) = rows2use(1) - 1;
	Nacceptor = rows2use(2)-rows2use(1);
    US = zeros(Nacceptor,size(U,2));
	for(i=1:Nacceptor) % explicit multiplication by columns to save mem. vs. explicit mat. mult of U*S
		US(i,:) = U(i+rows2use(1),:).*S; 
    end
    
%	clear S; % save on all the memory we can
%	Uprime = (U(cols2use(1):cols2use(2),:))';
%	clear U;

	% write covar reconstruction a chunk at a time

	Ndonor = 1+cols2use(2)-cols2use(1);
	Ndonor_pfacts = factor(Ndonor);
	chunk_size = prod(Ndonor_pfacts(1:ceil(length(Ndonor_pfacts)/2)));

    U = U(cols2use(1):cols2use(2),:); % trim U to be only what we need - hopefully will save memory
    
    if(Z_flag == 1) % write Z-transform
        Nd = size(U,2);
        var_A = sigma_A^2;
        var_B = sigma_B^2;
        nn = -Nd*var_A*var_B;
        
        Adiag = sum(US.^2,2);
        
    	for(i=1:chunk_size:Ndonor)
            USprime_chunk = (U(i:(i+chunk_size-1),:).*repmat(S,chunk_size,1))';
            Bdiag_chunk = sum(USprime_chunk.^2,1);
            rows2write = (US*USprime_chunk)./ ...
                sqrt(nn + var_B*repmat(Adiag,1,size(USprime_chunk,2)) + var_A*repmat(Bdiag_chunk,size(US,1),1)); 
            bytes_written = fwrite(pipe_outfile, rows2write, 'float32');
        end
    elseif(Z_flag == 2) % write mrs transform
        Amax = max(max(abs(US*V'),min(sigma_A,sigma_B)/10000),[],2);
        
        for(i=1:chunk_size:Ndonor)
            USprime_chunk = (U(i:(i+chunk_size-1),:).*repmat(S,chunk_size,1))';
            Bmax_chunk = (max(max(abs(USprime_chunk'*V'),min(sigma_A,sigma_B)/10000),[],2)').^(-1);
            rows2write = (US*USprime_chunk)./(exp(abs(log(Amax*Bmax_chunk))));
            bytes_written = fwrite(pipe_outfile, rows2write, 'float32');
        end
    else
        for(i=1:chunk_size:Ndonor)
            USprime_chunk = (U(i:(i+chunk_size-1),:).*repmat(S,chunk_size,1))';
            rows2write = US*USprime_chunk;
            bytes_written = fwrite(pipe_outfile, rows2write, 'float32');
        end
    end

else
	fwrite(pipe_outfile,C(:),'float32'); % have covar result itself in C, so just write it out
end
fclose(pipe_outfile);