Merge remote-tracking branch 'origin/master'

int_prob.m

@@ -5,28 +5,35 @@ function p = int_prob(l, T, x0, direc, ind_t, bp, ind_delta, lb, T_mov)
 % bp        ... boundary position at t==ind_t
 % lb        ... distance from boundary not to take into account
 
+ss = T{1}.system_size;
 delta_x0 = diff(x0);
 p = 0;
+
+if lb ~= 0 && abs(direc)==2
+    disp('Jump within same phase not implemented for lb!=0');
+    p = nan;
+    return;
+end
+
 for i = 1:length(delta_x0)
     x = (x0(i)+x0(i+1))/2;
     % 1. cond.: corr. starting point?  2. cond: jumped outside of domain?
     corr_starting_point = direc*x < direc*bp-lb;
     if abs(direc) == 1 % jump across the boundary?
         corr_end_point = direc*(x-l) > direc*(bp-T{i}.v*T{1}.t(ind_delta+1))+lb;
-    elseif abs(direc) == 2 % jump within the same phase?
-        if lb ~= 0
-            disp('Jump within same phase not implemented for lb!=0');
-            break;
+    elseif direc == 2 % jump left -> left
+        if l < 0 && x-l < bp; corr_end_point = 1; else; corr_end_point=0; end
+        if l > 0 % jump to left
+            corr_end_point = 1;
+            if x-l < 0; x = l-x; end % reflecting boundary if jump out left
         end
-        l = abs(l); % both directions need to be considered.
-        if (x < bp)
-            if (x + l < bp);    right = 1;  else; right = 0;    end
-            if (x - l > 0);   left = 1;   else; left = 0;     end
-        else
-            if (x + l < T{1}.system_size);  right =1;   else; right = 0;    end
-            if (x -l > bp);     left = 1;   else; left = 0;     end
+    elseif direc == -2 % jump right -> right
+        if l > 0 && x-l > bp; corr_end_point = 1; else; corr_end_point=0; end
+        if l < 0 % jump to right
+            corr_end_point = 1;
+            % Reflecting boundary if jumping out on right side of system
+            if x-l > ss; x = 2*ss+l-x; end
         end
-        corr_end_point = left || right; % consider, if left or right works
     end
     if corr_starting_point && corr_end_point
         if nargin==8
@@ -40,17 +47,13 @@ for i = 1:length(delta_x0)
                     T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0, 0)*p2;
         elseif nargin==9
             if abs(direc) == 1
-                p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_delta+1, :), x-l);
-            elseif abs(direc) == 2
-                x_l = x-l; % left position
-                x_r = x+l; % right position
-                if x-l < 0; x_temp = l-x; end % -> jump reflected if necess.
-                if x+l > T{1}.system_size; x_temp = 2*T{1}.system_size-x-l; end
-                p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_delta+1, :),...
-                                   x_l) * left + ...
-                           interp1(T{j}.x, T{j}.sol(ind_delta+1, :),...
-                                   x_r) * right;
+                x_ev = x-l;
+            elseif direc == 2 % left -> left
+                if l > 0 && x-l < 0; x_ev = l-x; else; x_ev = x-l; end
+            elseif direc == -2 % right -> right
+                if l < 0 && x+l > ss; x_ev = 2*ss+l-x; else; x_ev = x-l; end
             end
+            p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_delta+1, :), x_ev);
             p2 = (p_i(i)+p_i(i+1))/2;
             p_sol = @(j) interp1(T_mov.x, T_mov.sol(ind_t, :), x);
 %             p_sol2 = T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0, 0);

prob_laplace.m

@@ -6,6 +6,9 @@ e = @(chi) sqrt(3/8*(chi-2));
 t = linspace(0.0, 0.001, 11);
 % t=0.01;
 % t = logspace(-4, -2, 20);
+% Folders
+next = '/home/hubatsch/Nextcloud/';
+mac_next = '/Users/hubatsch/Nextcloud/';
 %% Comparison with Python
 T1 = Ternary_model(0, 'FRAP', {-0.1, b(7/3, 10^-12), 0.5, 0.4,...
                                      0, 1, 300, 7, 0, 'Constituent', 0},...
@@ -65,8 +68,8 @@ for i = 1:length(fac)
                               e, e_g, u_g, 300, 7, 0, 'Constituent'},...
                               t, 2);
 %     T_prec(i) = Ternary_model(2, 'FRAP', {-1, b(7.7/3, 10^-5), u0, ...
-%                               e(7.7/3), e_g, u_g, 300, 7, 0, 'Constituent'},...
-%                               t, 0.2);
+%                               e(7.7/3), e_g, u_g, 300, 7, 0,...
+%                               'Constituent'}, t, 0.2);
     T_prec(i).solve_tern_frap()
 end
 
@@ -86,9 +89,9 @@ Gi = 1;%(1-Tpt1.phi_t(1))/(1-Tpt1.phi_t(1));
 Go = Tpt1.phi_t(end)/Tpt1.phi_t(1)*2;
 %%
 for i = 1:1%length(prec)
-    Tga1(i) = Ternary_model(2, 'FRAP', {-1, b(7/3, 10^-12), Tpt1(i).phi_t(1),...
-                                        0, -0.83, Gi, 300, 7, 0,...
-                                        'Constituent'},...
+    Tga1(i) = Ternary_model(2, 'FRAP', {-1, b(7/3, 10^-12), ...
+                                        Tpt1(i).phi_t(1), 0, -0.83, Gi,...
+                                        300, 7, 0,'Constituent'},...
                             t, prec(i));
     Tga1(i).x = Tpt1(i).x;
     Tga1(i).solve_tern_frap();
@@ -230,7 +233,7 @@ sum((T{1}.sol(1, 1:end-1)+T{1}.sol(1, 2:end))/2.*diff(T{1}.x))
 % figure; hold on;
 plot(ls, p/N);
 %% Write to file
-cd /Users/hubatsch/Nextcloud/Langevin_vs_MeanField/Data_Figs_FokkPla/jump_length/
+cd([next_mac, 'Langevin_vs_MeanField/Data_Figs_FokkPla/jump_length/']);
 csvwrite('jump_length_7_7_lb01.csv', ls)
 csvwrite('prob_7_7_lb01.csv', p);
 
@@ -427,9 +430,9 @@ for j = 1:n_T
     end
     toc
 end
-
-% csvwrite('/home/hubatsch/Nextcloud/jump_length_mov_bound_in_out.csv', [-ls', sum(p, 2)])
-% csvwrite('/home/hubatsch/Nextcloud/jump_length_mov_bound_out_in.csv', [ls', sum(p, 2)])
+next = '/home/hubatsch/Nextcloud/';
+% csvwrite([next, 'jump_length_mov_bound_in_out.csv'], [-ls', sum(p, 2)])
+% csvwrite([next, 'jump_length_mov_bound_out_in.csv', [ls', sum(p, 2)])
 %% %%%%%%%%%%%%%%%%%%%%%% FRAP TIME COURSE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 t = linspace(0, 5, 300);
@@ -492,43 +495,55 @@ bp = 5;
 direc = -2; % 1: jump from left to right, -1: jump from right to left.
 
 %% 
-[ls, p_i_i] = frap_distro(t, params, 2, bp, 1);
-[~, p_o_o] = frap_distro(t, params, -2, bp, 1);
-[~, p_i_o] = frap_distro(t, params, 1, bp, 1);
-[~, p_o_i] = frap_distro(t, params, -1, bp, 1);
+[ls_i_i, p_i_i] = frap_distro(t, params, 2, bp, 1);
+[ls_o_o, p_o_o] = frap_distro(t, params, -2, bp, 1);
+[ls_i_o, p_i_o] = frap_distro(t, params, 1, bp, 1);
+[ls_o_i, p_o_i] = frap_distro(t, params, -1, bp, 1);
 %%
 p_i_i(isnan(p_i_i)) = 0;
 p_o_o(isnan(p_o_o)) = 0;
 p_i_o(isnan(p_i_o)) = 0;
 p_o_i(isnan(p_o_i)) = 0;
 %%
+i = 2;
+figure; hold on;
+plot(ls_i_i, p_i_i(:, i));
+plot(ls_o_o, p_o_o(:, i));
+plot(ls, p_i_o(:, i));
+plot(ls, p_o_i(:, i));
+%%
 figure; hold on;
-plot(ls, sum(p_i_i, 2));
-plot(ls, sum(p_o_o, 2));
-plot(ls, sum(p_i_o, 2));
-plot(ls, sum(p_o_i, 2));
-save('/data/biophys/hubatsch/MatlabWorkspaces/prob_laplace_190121.mat');
+plot(ls_i_i, sum(p_i_i, 2));
+plot(ls_o_o, sum(p_o_o, 2));
+plot(ls_i_o, sum(p_i_o, 2));
+plot(ls_o_i, sum(p_o_i, 2));
+% save('/data/biophys/hubatsch/MatlabWorkspaces/prob_laplace_190121_2.mat');
 %% Save time points for comparison with python
 csvwrite('lb_in_out.csv', ls)
 csvwrite('prob_in_out.csv', p);
 %%
-KL = zeros(1, 100);
-for i = 1:68
-i_t = 1+i;
-j_t = 2+i;
 integrand = @(p, q) abs(p.*log(p./q));
 integrate = @(x, y) nansum(diff(x)'.*(y(2:end)+y(1:end-1))/2);
-n_i = @(p1, p2, p3, p4) abs(integrate(ls, p1) + integrate(ls, p2)+...
-                            integrate(ls, p3) + integrate(ls, p4));
+KL = zeros(1, 50);
+for i = 1:42
+i_t = 1+i;
+j_t = 2+i;
+n_i = @(ii, oo, io, oi) abs(integrate(ls_i_i, ii) + integrate(ls_o_o, oo)+...
+                            integrate(fliplr(ls_i_o), flipud(io)) +...
+                            integrate(ls_o_i, oi));
 N_i = n_i(p_i_i(:, i_t), p_o_o(:, i_t), p_i_o(:, i_t), p_o_i(:, i_t));
 N_j = n_i(p_i_i(:, j_t), p_o_o(:, j_t), p_i_o(:, j_t), p_o_i(:, j_t));
-
-KL(i) = integrate(ls, integrand(p_i_i(:, i_t)/N_i, p_i_i(:, j_t)/N_j))+...
-     integrate(ls, integrand(p_o_o(:, i_t)/N_i, p_o_o(:, j_t)/N_j))+...
-     integrate(ls, integrand(p_i_o(:, i_t)/N_i, p_o_i(:, j_t)/N_j))+...
-     integrate(ls, integrand(p_o_i(:, i_t)/N_i, p_i_o(:, j_t)/N_j));
+i_i = integrate(ls_i_i, integrand(p_i_i(:, i_t)/N_i,...
+                flipud(p_i_i(:, j_t))/N_j));
+o_o = integrate(ls_o_o, integrand(p_o_o(:, i_t)/N_i,...
+                flipud(p_o_o(:, j_t))/N_j));
+i_o = integrate(fliplr(ls_i_o), integrand(flipud(p_i_o(:, i_t))/...
+                N_i, p_o_i(:, j_t)/N_j));
+o_i = integrate(ls_o_i, integrand(p_o_i(:, i_t)/N_i, p_i_o(:, j_t)/N_j));
+KL(i) = i_i + o_o + i_o + o_i;
 end
-plot(t(2:101), abs(KL)/0.02);
+
+plot(t(2:51), KL/0.02);
 % ylim([0, 1])
 %%
 figure; hold on;
@@ -539,7 +554,9 @@ plot(ls, p_o_i(:, j_t)/N_j, 'b');
 %%
 plot(T_mov.t, s_dot)
 hold on; xlim([0, 2]);
-plot(t(2:101), 0.1908/0.1885*0.61*abs(KL)/0.02);
+plot(t(2:101), 0.617*abs(KL)/0.02);
+%%
+plot(t(2:51), 0.557*abs(KL)/0.04);
 %% %%%%%%%%%%%%%%%%%%% INTEGRATION FUNCTION DEFINITIONS %%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -560,7 +577,10 @@ end
 T_mov = Ternary_model(0, 'FRAP', params, t, 0.2);
 T_mov.solve_tern_frap();
 ls = -abs(direc)/direc*(0.001:0.04:4);
-n_T=70;
+if abs(direc) == 2 % jumps within same phase
+    ls = sort([ls, -ls]); % need pos. and neg. l, since direction not set
+end
+n_T=45;
 p = nan(length(ls), n_T);
 for j = 1:n_T
     tic
@@ -597,6 +617,6 @@ dudx = diff(Temp.sol(t_ind, :))./diff(Temp.x);
 x = (Temp.x(1:end-1)+Temp.x(2:end))/2;
 sol = (Temp.sol(t_ind, 1:end-1)+Temp.sol(t_ind, 2:end))/2;
 
-[~, flux ,~, g0, pt] = flory_hugg_pde(x, 0, sol, dudx, Temp.a, Temp.b, Temp.e, ...
-               Temp.u0, Temp.e_g0, Temp.u_g0, Temp.v, Temp.mode);
+[~, flux ,~, g0, pt] = flory_hugg_pde(x, 0, sol, dudx, Temp.a, Temp.b, ...
+                Temp.e, Temp.u0, Temp.e_g0, Temp.u_g0, Temp.v, Temp.mode);
 end
\ No newline at end of file