clear all;
close all;
clc;
%
COST = importdata('Outputs/cost.txt');
cost=COST(:,8); 
[min2,row2] = min(COST(:,8));
%
q0 =COST(row2,1)/90*(pi/2)  ;
q1 =COST(row2,2)/90*(pi/2)   ;
q2 =COST(row2,3)/360*(2*pi);
%
% --------------------3D rotation illustration------------------------------      
b = [0; 0; -1];                                                          % beam direction
n0 = [0; sin(q0); cos(q0)];                                      % initial surface tilt-to direction
a = [sin(q1)*cos(q2); sin(q1)*sin(q2); cos(q1)];     % rotation axis

figure(1);                                                                 % 3D rotation illustration of the ion momentum mixing 
hold on;
plot3([0 -b(1)], [0 -b(2)], [0 -b(3)], 'b','LineWidth',5); 
plot3([0 n0(1)], [0 n0(2)], [0 n0(3)], 'g','LineWidth',5);
plot3([0 a(1)], [0 a(2)], [0 a(3)], 'r','LineWidth',5);
legend('Ion beam', 'Surface normal', 'Rotation axis','Location','EastOutside')

mesh = 45; dphi = 2*pi/mesh;                              

for phi = 0 : dphi : 2*pi-dphi;
      new     =   [cos(q1)*sin(q1)*cos(q2);cos(q1)*sin(q1)*sin(q2);cos(q1)^2]+...
                        sin(q1)*(cos(phi)*[-cos(q1)*cos(q2);-cos(q1)*sin(q2);sin(q1)]...
                        + sin(phi)*[sin(q2);-cos(q2);0]);

       cosine = sin(q0)*sin(q1)*[cos(q1)*sin(q2)*(1-cos(phi))-cos(q2)*...
                       sin(phi)]+cos(q0)*[1-sin(q1)^2*(1-cos(phi))];

       if (cosine > 0)
                 plot3([0 new(1)], [0 new(2)], [0 new(3)], 'k');
       else
                 plot3([0 new(1)], [0 new(2)], [0 new(3)], 'k--');
       end
end
title(sprintf('\\theta_0=%g \\theta_1=%g \\theta_2=%g',q0,q1,q2));
set(gca, 'FontSize', 16,'FontWeight','bold');
view(70,30);
axis equal; axis vis3d;         
set(gca,'xTick',[0:0.5:1]);
set(gca,'yTick',[-0.5:0.5:0.5]);
set(gca,'zTick',[0:0.5:1]);
hold off;
%print -dtiff -r400 Outputs/3D-rotation

%% --------------Directional cosine in one-linear forms ------------------------

figure(2)
%hold on;
mesh=180; dphi=2*pi/mesh;                             % Mesh grid defined
phi=0:dphi:2*pi-dphi;
cosine = sin(q0)*sin(q1)*[cos(q1)*sin(q2)*(1-cos(phi))-cos(q2)*...
                sin(phi)]+cos(q0)*[1-sin(q1)^2*(1-cos(phi))];
plot(phi,cosine,'k','LineWidth',3)
%fitting with directional cosine to a diurnal cycle, 
%to obtain parameters: A B C
xdata=phi; ydata=cosine;
fun = @(a,xdata)a(1)*cos(xdata+a(2))+a(3);
x0 = [0.1, 2, 0.8];
lb = [0.0, 0.1, -2 ];
ub = [5,   10,  15];
xfit = lsqcurvefit(fun,x0,xdata,ydata,lb,ub)
%
phi_f= linspace(xdata(1),xdata(end));
cosine_f=fun(xfit, phi_f);
%plot(phi_f,cosine_f,'g-','LineWidth',2);

xlabel('\phi = {\omega}\itt'); ylabel('\itf(\itt) {\equiv} cos{\theta}(\itt)');
ylim([0 1]);
xlim([0 2*pi]);

set(gca,'xTick',[0:1:2*pi]);
set(gca,'yTick',[0:0.5:1]);
set(gca, 'FontSize', 20,'FontWeight','bold')
set(gcf,'position',[0,0,500,300])
width=500;%¿í¶È,ÏñËØÊý 
height=300

grid on; box on;
hold off;
%print -dtiff -r400 Outputs/f-factor
%
fid_Cosine= fopen('Outputs/Cosine_Params.txt','w')                    %Output the DOFs & parametrs of the fitted sun circle to a file.  
fprintf(fid_Cosine,'%s','The optimized DOFs of SSS are: ')
fprintf(fid_Cosine,'%30.15f',q0,q1,q2);fprintf(fid_Cosine,'\n');
fprintf(fid_Cosine,'%30.15f',COST(row2,1),COST(row2,2),COST(row2,3))
fprintf(fid_Cosine,'\n')
fprintf(fid_Cosine,'%s','The fitted Parameters are: ')
fprintf(fid_Cosine,'%30.15f',xfit)
fprintf(fid_Cosine,'\n')
fprintf(fid_Cosine,'\n')
fclose(fid_Cosine);  

%% ------- Average damage and excess designed by SSS-----------------------

% Paremeters inputs for SSS
IP = fopen('input.txt','r');                                    % scan all data extracting variables
instring = fgets(IP);
E0 = sscanf(instring(1:10), '%e')                         % E0, monochromatic kinetic energy of ion (MeV) 
instring = fgets(IP);
Zmin = sscanf(instring(1:10), '%e')                     % Zmin, minimum boundary of desired depth region (¦Ìm)         
instring = fgets(IP);
Zmax = sscanf(instring(1:10), '%e')                    % Zmax, Maximum boundary of desired depth region (¦Ìm)

E0; Zmin; Zmax;                                                   
E = importdata('Outputs/total-excess.txt');
D= importdata('Outputs/total-damage.txt');
z = (0.015 : 0.005 : 2.10) ;                                     %%% total depth 

minZ=find(abs(z-Zmin)<=1e-4);
maxZ=find(abs(z-Zmax)<=1e-4);
Zdel = 0.005;
Z = (Zmin : Zdel : Zmax);                                   

HH = 0;
sumHH =0;
sumDD =0;
sumDDZ =0;
sumEEZ=0;
aveDD = 0;
aveDDZ = 0;
aveEE =0;
aveEEZ=0;
AA = 0;

figure(3);                                                              %% Average damage plot
set(gcf,'position',[0,0,700,350])
yyaxis left;
 for nn=1:length(phi)
            theta(nn) = acos(cosine(nn))*180/pi;
            theta(nn)=roundn(theta(nn),0)
           
            if((theta(nn))<=89 )            
                    HH = heaviside(cosine(nn)) ;  
                    AA = AA + dphi* HH.* cosine(nn);
                   size(z)
                   size(E(:,theta(nn)+1))
                   %  plot(z, E(:,theta(nn)+1)'); hold on;
                   if(mod(nn,2) == 0)
                      plot(z(2:end-1), D(2:end-1,theta(nn)+1)','-','color',[0.6 0.6 1],'LineWidth',0.05); hold on;
                   end
                   sumDD    = sumDD   + dphi*HH.*cosine(nn)* D(    :   ,          theta(nn)+1)' ; 
                   sumDDZ  = sumDDZ + dphi*HH.*cosine(nn)* D(minZ:maxZ,theta(nn)+1)' ; 
                   sumHH    = sumHH   + dphi*HH.*cosine(nn)* E(     :   ,          theta(nn)+1)' ; 
                   sumEEZ    = sumEEZ + dphi*HH.*cosine(nn)* E(minZ:maxZ,theta(nn)+1)' ;
            end 
            aveDD   = sumDD/AA;
            aveDDZ = sumDDZ/AA;
            aveEE     = sumHH/AA;
            aveEEZ   = sumEEZ/AA;
end     
plot([0 2], [8 8], 'k--','LineWidth',1.2);
plot([Zmin Zmin],[0 40],'k--','LineWidth',1.2);
plot([Zmax Zmax],[0 25],'k--','LineWidth',1.2);
plot(z(2:end-1), aveDD(2:end-1),'b','LineWidth',2.5);
plot(Z,aveDDZ,'b','LineWidth',5);
plot(z,D(:,1),'b-.','Linewidth',2);
 
axis([0 2 -0.2 40]);
set(gca,'ytick',[-0.0:10:40]);
set(gca,'xtick',[0:0.5:2.0]);
xlabel('Depth (¦Ìm)');
ylabel('Vacancy / ion / nm' ,'FontSize', 14,'FontWeight','bold');
set(gca,'YColor','b');
set(gca,'Position',[0.12,0.20,0.77,0.65])
set(gca, 'FontSize', 16,'FontWeight','bold')
box on;
yyaxis right;
   set(gca,'yticklabel',[])
   set(gca,'YColor','k');
hold off;
%print -dtiff -r400 Outputs/average-Dam
%%
figure(4);                                                             % Average excess plot
hold on;
set(gcf,'position',[0,0,700,400]);
 for nn=1:1:length(phi)
       theta(nn) = acos(cosine(nn))*180/pi;
       theta(nn)=roundn(theta(nn),0);
       if((theta(nn))<=89 )            
               if(mod(nn,2) == 0)                              % not all the channels are plotted
                      plot(z(2:end-1), E(2:end-1,theta(nn)+1)','-','color',[1 0.6 0.6],'LineWidth',0.05); 
               end
       end 
 end 
plot([0 2], [0 0], 'k--','LineWidth',1.2);
plot([Zmin Zmin],[-0.003 0.002],'k--','LineWidth',1.2);
plot([Zmax Zmax],[-0.001 0.002],'k--','LineWidth',1.2);
plot(z,E(:,1),'r-.','Linewidth',2);
plot(z(2:end-1), aveEE(2:end-1),'k','LineWidth',2.5);
plot(Z,aveEEZ,'k','LineWidth',5);

%set(gca,'YColor','k');
axis([0 2 -0.003 0.002]);
set(gca,'ytick',[-0.003:0.001:0.002]);
set(gca,'xtick',[0:0.5:2.0]);
xlabel('Depth (¦Ìm)');
ylabel('I-V excess / ion / nm' ,'FontSize', 14,'FontWeight','bold');
set(gca,'Position',[0.12,0.20,0.77,0.65])
set(gca, 'FontSize', 16,'FontWeight','bold')
box on;
hold off;
%print -dtiff -r400 Outputs/average-Exc
%% RMSE calculation

fid_RMSR= fopen('Outputs/RMSE.txt','w')  
% RMSE of damage, the reference value is the average damage in the area of interest
Rp=1.8;                                                % depth of ion hit contours estimated from I-V excess concentration
pos=find(abs(z-Rp)<=1e-4) ;              
muD= mean(D(1:pos,1));                     % ref: average damage in ion hit contours for normal irradiation
lenD= pos;
%
RMSE_D=  sqrt(sum( (D(1:pos,1)-muD).^2)./lenD);    % RMSE of Damage for normal irradiation
         
muDZ= mean(aveDDZ) ;                    % ref: average damage (DDZ) in desired depth region under SSS
%
lenDZ=length(aveDDZ);
RMSE_DZ=  sqrt(sum( (aveDDZ-muDZ).^2)./lenDZ);    % RMSE of Damage under SSS

fprintf(fid_RMSR,'\n');
fprintf(fid_RMSR,'%s','The RMSE of DAMAGE before and after using SSS are:');
fprintf(fid_RMSR,'%10.2f',RMSE_D,RMSE_DZ);
fprintf(fid_RMSR,'\n');

% RMSE of I-V excess, the reference value is ZERO
muE= 0;                                                                    % reference value: 0 
lenE= pos;
RMSE_E=  sqrt(sum( (E(1:pos,1)-muE).^2)./lenE);    % RMSE of Damage for normal irradiation
         
muHZ= 0 ;                                                                %  reference value:  0 
lenHZ=length(aveEEZ);
RMSE_EZ=  sqrt(sum( (aveEEZ-muE).^2)./lenHZ);    % RMSE of Damage under SSS

fprintf(fid_RMSR,'%s','The RMSE of EXCESS before and after using SSS are:');
fprintf(fid_RMSR,'%10.2e',RMSE_E, RMSE_EZ);
fprintf(fid_RMSR,'\n');
fclose(fid_RMSR); 
%
fprintf('%s','Displaying is done')