% porcelain_rga_data_analaysis.m
% Read asci data from RGA scan
% 2010-06-15, D. Coyne
dir='C:\L2\Sys\Uhv\material qualification\porcelain\';
outfile='porcelainized SS outgassing.txt';

% Material
% 29 plates of 18 gauge stainless steel, each 15 inches x 30 inches
% porcelainized on both sides with Ferro Corp. PE RTU GC Gloss Blk RM108 frit
area=29*15*30*2*2.54^2;

% Read
filename='F51210.txt';
fid=fopen(char(strcat(dir,filename)),'r');
FILE_DATA = (fread(fid,'char=>char')).';
ind = strfind(FILE_DATA,sprintf('\n'));

% echo header
ofid=fopen(char(strcat(dir,outfile)),'w');
fprintf(ofid,'%s',FILE_DATA(1:ind(22)));

% read RGA data (AMU,Amps)
[data,count]=sscanf(FILE_DATA(ind(22)+1:end),'%f %*s %f %*s');
st = fclose(fid);
n=count/2;
d=reshape(data,2,n);
amu=d(1,:);
amp=d(2,:);
%
% plot
figure(1);
hold off;
semilogy(amu,amp,'g-');
title('Porcelainized Stainless Steel (Ferro Corp. PE RTU GC Gloss Blk RM108 frit)');
axis([0 100 1e-17 1e-11]);
xlabel('AMU');
ylabel('amp');
grid on;
hold on;

% Calibrated Leak
LR=2.36e-10; %torr-liters/s
% Calibration factor
calAmp=1.8e-14;
cal=LR/calAmp; %torr-liter/sec/Amp

% AMU Shift
% Find the AMU associated with the Ar peak at 44
ind1=find(amu > 43);
ind2=find(amu > 45);
[a,b]=max(amp(ind1(1):ind2(1)));
shift=amu(ind1(b))-44;

% Hydrocarbon Outgassing Rate
ind41=find(amu >= (41+shift));
ind43=find(amu >= (43+shift));
ind53=find(amu >= (53+shift));
ind55=find(amu >= (55+shift));
ind57=find(amu >= (57+shift));
amp41=amp(ind41(1));
amp43=amp(ind43(1));
amp53=amp(ind53(1));
amp55=amp(ind55(1));
amp57=amp(ind57(1));
HC=(amp41+amp43+amp53+amp55+amp57)*cal;
ind80=find(amu >= 80+shift);
ampNoiseLevel=mean(amp(ind80:end));
noiselimitHC=5*ampNoiseLevel*cal;

% find spectrum peaks
[pks1,locs1]=findpeaks(amp,'minpeakdistance',3);
[pks2,locs2]=findpeaks(amp,'minpeakheight',5*ampNoiseLevel);
[pks3,locs3]=findpeaks(amp,'threshold',ampNoiseLevel/10);
locs=intersect(intersect(locs1,locs2),locs3);
amuShiftErrors=ceil(amu(locs))-amu(locs);

% plot peaks on RGA scan
semilogy(ceil(amu(locs)),amp(locs),'ko');

% linear fit to peak shift vs AMU
p=polyfit(amu(locs),amuShiftErrors,2);
amuShiftFit=p(1)*amu(locs).^2+p(2)*amu(locs)+p(3);
%amuShiftFit=p(1)*amu(locs)+p(2);
figure;
plot(amu(locs),amuShiftErrors,amu(locs),amuShiftFit);
title('Porcelain');
xlabel('AMU');
ylabel('AMU shift error');
grid on;

% Find spectrum values at key AMUs
gases=['H2 ';'He ';'CO ';'N2 ';'O2 ';'Ar ';'CO2';'H2O';'HC '];
H2=[1,2,3];
He=[2,4];
CO=[12,14,16,28,29,30];
N2=[14,28,29,30];
O2=[8,16,17,18,32,33,34,36];
Ar=[18,19,20,36,38,40];
CO2=[12,13,16,18,22,23,28,29,30,44,45,46,47];
H2O=[1,16,17,18,19];
a=union(union(union(union(union(union(union(H2,He),CO),N2),O2),Ar),CO2),H2O);
aShifted=a-p(1)*a.^2-p(2)*a-p(3);
n=length(a);
index=zeros(1,n);
for i=1:n
    tU=find(amu>=aShifted(i));
    t1=tU(1);
    tL=find(amu<=aShifted(i));
    if isempty(tL),
        index(i)=t1;
    else
        t2=tL(1);
        if abs(amu(t1)-aShifted(i))<abs(amu(t2)-aShifted(i)),
           index(i)=t1;
        else
            index(i)=t2;
        end
    end
end
b=amp(index);

index=zeros(1,n);
nsteps=2;
for i=1:n
    tU=find(amu>=aShifted(i));
    index(i)=tU(1);
    amax=0;
    aind=index(i);
    for j=-nsteps:nsteps
        if (index(i)+j)>0 & amp(index(i)+j)>amax,
            aind=index(i)+j;
            amax=amp(index(i)+j);
        end
    end
    index(i)=aind;
end
b=amp(index);

% plot key AMU values on RGA scan
figure(1);
semilogy(a,b,'rx');

% Cracking pattern matrix
% AMUS = 1,2,3,4,8,12,13,14,16,17,18,19,20,22,23,28,29,30,32,33,34,36,38,40,44,45,46,47
% Gases = H2, He, CO, N2, O2, Ar, CO2, H2O
nAMU=28;
nGAS=8;
c=zeros(nAMU,nGAS);
c(1,1)=.027;
c(2,1)=1;
c(3,1)=.0031;
c(2,2)=.0012;
c(4,2)=1;
c(5,3)=.000001;
c(6,3)=.035;
c(8,3)=.014;
c(9,3)=.014;
c(16,3)=1;
c(17,3)=.012;
c(18,3)=.002;
c(8,4)=.09;
c(16,4)=1;
c(17,4)=.0071;
c(18,4)=.000014;
c(5,5)=.000013;
c(9,5)=.14;
c(10,5)=.000052;
c(11,5)=.00028;
c(19,5)=1;
c(20,5)=.00074;
c(21,5)=.0038;
c(22,5)=.000023;
c(11,6)=.00071;
c(12,6)=.00016;
c(13,6)=.05;
c(22,6)=.0036;
c(23,6)=.00068;
c(24,6)=1;
c(4,7)=.000005;
c(5,7)=.000005;
c(6,7)=.063;
c(7,7)=.00063;
c(9,7)=.16;
c(11,7)=.000088;
c(14,7)=.0052;
c(15,7)=.000012;
c(16,7)=.15;
c(17,7)=.0015;
c(18,7)=.00029;
c(25,7)=1;
c(26,7)=.012;
c(27,7)=.0038;
c(28,7)=.000034;
c(1,8)=.001;
c(9,8)=.0307;
c(10,8)=.2701;
c(11,8)=1;
c(12,8)=.0019;

% Total Ionization Cross Sections
% Gases = H2, He, CO, N2, O2, Ar, CO2, H2O
cs=[0.42,0.14,1.07,1,1.02,1.19,1.36,1]';
for i=1:nGAS
    c(:,i)=c(:,i).*cs(i);
end

% Solve for gas outgassing rates
p=cal*pinv(c)*b'; % torr-liter/sec
p=vertcat(p,HC)/area;
fprintf(ofid,'Outgassing rates by gas species (torr-liter/sec/cm^2):\n');
for i=1:nGAS+1
    fprintf(ofid,'%s %e \n',gases(i,:),p(i));
end