TTB_rotate.m 2.14 KB
 Jürgen Matzka committed Apr 27, 2018 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 ``````function [p]=TTB_rotate(date,PATH) % Read data TTB1=read_data(PATH,'TTB1',date); TTB0=read_data(PATH,'TTB0',date); % Calculate F from rotated values %F=@(V,p) sqrt((V.HNvar.*Sc(V.T1,p(3)).*cosd(p(1))).^2+(V.HEvar.*Sc(V.T1,p(4))-V.HNvar.*Sc(V.T1,p(3)).*sind(p(1))).^2+(V.Zvar.*Sc(V.T1,p(5))).^2)+p(2); %F=@(V,p) sqrt((V.HNvar.*cosd(p(1))).^2+(V.HEvar-V.HNvar.*sind(p(1))).^2+(V.Zvar.^2))+p(2); %F=@(V,p) sqrt(((V.HNvar-V.HEvar./sind(p(1))).*cosd(p(1))).^2+(V.HEvar).^2+(V.Zvar.^2)); F=@(V,p) sqrt( ((V.HNvar-V.HEvar.*sind(p(1)))./cosd(p(1))).^2 + (V.HEvar).^2 + (V.Zvar.^2) ) +p(2); % Function for misfit value O=@(TTB1,p,n) nansum(abs(F(TTB1,p)-TTB1.Fsc).^n)^(1/n); % Fit n=2; % Least-squares p=fminsearch(@(p) O(TTB1,p,n),[0 0]); O(TTB1,p,n) % Plot clf subplot(4,1,1) plot(TTB1.time,TTB1.Fsc-TTB1.Fcalc,'.') title('Original \Delta F') subplot(4,1,2) plot(TTB1.time,TTB1.Fsc-F(TTB1,p),'.') title('Rotated \Delta F') subplot(4,1,3) plot(TTB1.time,-TTB1.T1) % Now temperature % Temperature scaling Sc=@(V,T,p) V.*(p(1).*T.^2+ones); %Sc=@(V,T,p) V+p(1).*T; %Sc=@(V,T,p) V; % Calculate F from rotated values r=p(1:2); F=@(V,p) sqrt( ((Sc(V.HNvar,V.T1,p(1))-Sc(V.HEvar,V.T1,p(1)).*sind(p(1)))./cosd(p(1))).^2 + (Sc(V.HEvar,V.T1,p(1))).^2 + (Sc(V.Zvar,V.T1,p(1)).^2) ); %F=@(V,p) sqrt((Sc(V.HNvar,V.T1,p(1)).*cosd(r(1))).^2+(Sc(V.HEvar,V.T1,p(1))-Sc(V.HNvar,V.T1,p(2)).*sind(r(1))).^2+(Sc(V.Zvar,V.T1,p(3))).^2)+r(2); %F=@(V,p) sqrt((Sc(V.HNvar,V.T1,p(1)).*cosd(r(1))).^2+(Sc(V.HEvar,V.T1,p(1))-Sc(V.HNvar,V.T1,p(1)).*sind(r(1))).^2+(Sc(V.Zvar,V.T1,p(1))).^2)+r(2); % Function for misfit value O=@(TTB1,p,n) nansum(abs(F(TTB1,p)-TTB1.Fsc).^n)^(1/n); % Fit n=2; % Least-squares p=fminsearch(@(p) O(TTB1,p,n),[0 0 0]); O(TTB1,p,n) % Plot subplot(4,1,4) plot(TTB1.time,TTB1.Fsc-F(TTB1,p),'.') title('Rotated \Delta F') figure plot(TTB1.time,TTB0.HNvar-(Sc(TTB1.HNvar,TTB1.T1,p(1)).*cosd(r(1)))); p=[r p]; end``````