%% collect all sessions variables
clear 
Sessions={ ...
'Y:\Analysed Data\Temp_for_RF\OM201\Multi_ContinousWave' ...
'Y:\Analysed Data\Temp_for_RF\OM213\MultiContiWave' ...
'Y:\Analysed Data\Temp_for_RF\OM214\MultiContWave' ...
'Y:\Analysed Data\Temp_for_RF\OM805\CW_50s200s\CW_11142021' ...
'Y:\Analysed Data\Temp_for_RF\OM806\CW50s200s' ...
};
Sessions_animals=["OM201","OM213","OM214","OM805","OM806",];
fig=figure(); hold on;

allData=[];
for sesi=1:length(Sessions)
    load([Sessions{sesi} '\allVars.mat']);
    allData=[allData allVars];
    if ~isempty(allVars)
        scatter([allVars.power_W],[allVars.T_rise_1sec],'filled','MarkerFaceAlpha',0.7)
    end
    clearvars -except Sessions sesi allData fig
end
save(['allData.mat'],'allData','-v7.3')

powerlevels_W=[sort(unique([allData.power_W]))];
meanAll_T_rise_1sec=[];
for i=1:length(powerlevels_W)
    idx=find([allData.power_W]==powerlevels_W(i));
    xx=[allData.T_rise_1sec];
    meanAll_T_rise_1sec(i)=mean(xx(idx));   
end

plot(powerlevels_W,meanAll_T_rise_1sec,'-xk','linewidth',2)
title(['1sTrise@37.52W= ' num2str(meanAll_T_rise_1sec(find(powerlevels_W==37.52))) ' C/s']);
xlabel('Radiate Power (W)')
ylabel('1sec \Delta T (C)')
set(gca,'fontweight','bold','fontsize',12,'linewidth',1)
set(fig,'Renderer', 'painters')
saveas(gca,fullfile(cd,['1sTriseVSpower.pdf']))
saveas(gca,fullfile(cd,['1sTriseVSpower.svg']))
saveas(gca,fullfile(cd,['1sTriseVSpower.fig']))
close all


%% Calcualate SAR and EF
%SAR
fig=figure(); hold on;
HeatCapa=3540.20; %J/kg·°C
for sesi=1:length(Sessions)
    load([Sessions{sesi} '\allVars.mat']);
    if ~isempty(allVars)
        scatter([allVars.power_W],HeatCapa*[allVars.T_rise_1sec],'filled','MarkerFaceAlpha',0.7)
    end
    clearvars -except Sessions sesi allData fig HeatCapa
end

powerlevels_W=[sort(unique([allData.power_W]))];
meanAll_T_rise_1sec=[];
for i=1:length(powerlevels_W)
    idx=find([allData.power_W]==powerlevels_W(i));
    xx=[allData.T_rise_1sec];
    meanAll_T_rise_1sec(i)=mean(xx(idx));   
end

meanAll_SAR=HeatCapa*meanAll_T_rise_1sec;

plot(powerlevels_W,meanAll_SAR,'-xk','linewidth',2)
xlabel('Radiate Power (W)')
ylabel('Calculated SAR (W/Kg)')
set(gca,'fontweight','bold','fontsize',12,'linewidth',1)
set(fig,'Renderer', 'painters')
title(['SAR(37.52W)= ' num2str(meanAll_SAR(find(powerlevels_W==37.52))) ' W/Kg']);

saveas(gca,fullfile(cd,['SARValue.svg']))
saveas(gca,fullfile(cd,['SARValue.pdf']))
saveas(gca,fullfile(cd,['SARValue.fig']))

%EF
fig=figure(); hold on;
density=1043; %kg/m3
ElecCondutivity=1.12;  % S/m
for sesi=1:length(Sessions)
    load([Sessions{sesi} '\allVars.mat']);
    if ~isempty(allVars)
        scatter([allVars.power_W],sqrt(2*density*HeatCapa*[allVars.T_rise_1sec]/ElecCondutivity),'filled','MarkerFaceAlpha',0.7)
    end
    clearvars -except Sessions sesi allData fig HeatCapa density ElecCondutivity meanAll_SAR
end
powerlevels_W=[sort(unique([allData.power_W]))];
meanAll_EF=sqrt(2*density*meanAll_SAR/ElecCondutivity);
plot(powerlevels_W,meanAll_EF,'-xk','linewidth',2)

xlabel('Radiate Power (W)')
ylabel('Estimated EF (V/m)')
set(gca,'fontweight','bold','fontsize',12,'linewidth',1)
set(fig,'Renderer', 'painters')
title(['EF(37.52W)= ' num2str(meanAll_EF(find(powerlevels_W==37.52))) ' V/m'])

saveas(gca,fullfile(cd,['EF_estimate.svg']))
saveas(gca,fullfile(cd,['EF_estimate.pdf']))
saveas(gca,fullfile(cd,['EF_estimate.fig']))


save(['VariablesAt' datestr(now,'mmddyyHHMMSS') '.mat'],'-v7.3')

close all