function plotOptomech();
% plot TM to TM optomachanical TFs
%H=importdata('./suspensionTF/oplev_120314/120314_TypeA_TMp.txt');
H=importdata('./suspensionTF/typeA_20150221/typeA_TMp.txt');
btypeA_tf_pitch=squeeze(H(:,2).*exp(i.*H(:,3)/180*pi));
H=importdata('./suspensionTF/typeA_20150221/typeA_TMy.txt');
btypeA_tf_yaw=squeeze(H(:,2).*exp(i.*H(:,3)/180*pi));
freq=H(:,1);

p=parambLCGT('BRSE','positive','pitch');
c=299792458;
%p.ITMXChr=1/14e3;
%p.ETMXChr=1/7.5e3;
gITM=1-p.Larm*p.ITMXChr;
gETM=1-p.Larm*p.ETMXChr;
r=2/((gITM-gETM)+sqrt((gITM-gETM)^2+4))
Pintracav=380e3;  %intra-cavity power
k0=2*Pintracav*p.Larm/(c*(gETM*gITM-1));
kSOFT=k0*(gETM+1/r)
kHARD=k0*(gITM-1/r)

plotoptomechTF(freq,btypeA_tf_pitch,kSOFT,kHARD,'pitch',['Pincav=', num2str(Pintracav/1e3),'kW, gITM=',num2str(gITM,'%.2f'),', gETM=',num2str(gETM,'%.2f')]);
plotoptomechTF(freq,btypeA_tf_yaw,kSOFT,kHARD,'yaw',['Pincav=', num2str(Pintracav/1e3),'kW, gITM=',num2str(gITM,'%.2f'),', gETM=',num2str(gETM,'%.2f')]);

function plotoptomechTF(freq,H,kSOFT,kHARD,dof,tit);
I=1/(abs(H(end))*(2*pi*freq(end))^2);
fres=sqrt(1/(I*abs(H(1))))/(2*pi);
fit=filterdesign(freq,[],[fres 1e4],1/I);
k=(2*pi*fres)^2*I;
fSOFT=sqrt((k+kSOFT)/I)/(2*pi);
fHARD=sqrt((k+kHARD)/I)/(2*pi);

figure();
plotbode(freq,[optomech(H,0),optomech(H,kSOFT),optomech(H,kHARD),fit]);
legend('Type-A no radiation','Type-A SOFT','Type-A HARD','fitted');
title(['(I=', num2str(I,'%.2f'),'kgm^2, fmech=', num2str(fres,'%.2f'),'Hz, fSOFT=', num2str(fSOFT,'%.2f'),'Hz, fHARD=',num2str(fHARD,'%.2f'),'Hz)'])
xlim([1e-3,1e2]);
subplot(2,1,1);
ylabel('Abs [rad/(Nm)]');
ylim([1e-5,1e5]);
xlim([1e-3,1e2]);
set(gca,'YTick',logspace(-5,5,11));
title(['Torque on TM to angle of TM (bLCGT, ',dof,', ',tit,')']);

saveas(gcf,['./results/optomech/bLCGT_TM',dof,'_',tit,'.png'])
disp(['bLCGT, I=', num2str(I,'%.2f'),'kgm^2 ',tit]);
disp([dof,': fmech=', num2str(fres,'%.2f'),'Hz fSOFT=', num2str(fSOFT,'%.2f'),'Hz fHARD=',num2str(fHARD,'%.2f'),'Hz']);

function tf=optomech(H,k);
tf=H./(1+k*H);