TTB_rotate.m 2.14 KB
Newer Older
Jürgen Matzka's avatar
Jürgen Matzka committed
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