%%  Modified_London Penetration Depth
%%
%%  Created by Pei Ting-Hang on 15/5/2018.
%%  Copyright 2018. All rights reserved.

clear all;                           
                                     
T_Tc_x02=241;                        
T_Tc_y10=799;                        
T_Tc_x04=523;                        
T_Tc_y12=645;                        
T_Tc_x06=805;                        
T_Tc_y14=491;                        
T_Tc_x08=1087;                       
T_Tc_y16=337;                        
T_Tc_y18=183;                        
                                     
T_Tc_x=[279 301 318 339 361 386 416 448 487 529 554 583 613 630 643 663 704 723 738 773 794 810 833 848 872 897 917 933 956 994 1032 1078 1137 1179 1227];
T_Tc_y=[797 795 794 793 796 792 791 789 785 783 782 775 766 765 752 756 738 729 712 707 694 685 672 669 654 634 624 612 600 564 530 473 376 273 97];      
                                     
T_Tc_inv_x00=429;                    
T_Tc_inv_x02=709;                    
T_Tc_inv_x04=990;                    
T_Tc_inv_x06=1273;                   
T_Tc_inv_x08=1553;                   
T_Tc_inv_x10=1834;                   
T_Tc_inv_y00=1300;                   
T_Tc_inv_y02=1073;                   
T_Tc_inv_y04=846;                    
                                     
T_Tc_inv_x=[1732 1756 1790];         
T_Tc_inv_y=[1065 1114 1182];         
                                     
x_01=(T_Tc_x08-T_Tc_x06)/2.0         
y_01=(T_Tc_y12-T_Tc_y16)/4.0         
                                     
x_inv_01=(T_Tc_inv_x10-T_Tc_inv_x06)/4.0
y_inv_01=(T_Tc_inv_y00-T_Tc_inv_y04)/4.0
                                     
for i=1:length(T_Tc_x)               
    T_Tc_xreal(i)=(T_Tc_x(i)-T_Tc_x04)/x_01*0.10+0.40
    T_Tc_yreal(i)=(T_Tc_y12-T_Tc_y(i))/y_01*0.10+1.20
end                                  
vector_size=length(T_Tc_x);          
                                     
for i=1:length(T_Tc_inv_x)           
    T_Tc_xreal(vector_size+i)=(T_Tc_inv_x(i)-T_Tc_inv_x00)/x_inv_01*0.10
    T_Tc_yreal(vector_size+i)=sqrt(1.0/((T_Tc_inv_y00-T_Tc_inv_y(i))/y_inv_01*0.10))
end                                  
                                     
for i=1:length(T_Tc_xreal)           
    T_Tc_xreal10(i)=power(T_Tc_xreal(i),10.0)
    T_Tc_xreal8(i)=power(T_Tc_xreal(i),8.0)
    T_Tc_xreal6(i)=power(T_Tc_xreal(i),6.0)
    T_Tc_xreal5(i)=power(T_Tc_xreal(i),5.0)
    T_Tc_xreal4(i)=power(T_Tc_xreal(i),4.0)
    T_Tc_xreal3(i)=power(T_Tc_xreal(i),3.0)
    T_Tc_xreal2(i)=power(T_Tc_xreal(i),2.0)
    T_Tc_xreal_inv(i)=1.0/T_Tc_xreal(i)    
    T_Tc_yreal_inv(i)=1.0/T_Tc_yreal(i)    
    T_Tc_yreal_inv2(i)=power(1.0/T_Tc_yreal(i),2.0)
                                           
    delta_London_lambda2(i)=power(T_Tc_yreal(i),2.0)-(1.0/(1.0-power(T_Tc_xreal(i),4.0)))
    delta_London_lambda(i)=T_Tc_yreal(i)-sqrt(1.0/(1.0-power(T_Tc_xreal(i),4.0)))
    delta_T_Tc(i)=(1.0-power(T_Tc_xreal(i),4.0))-power(1.0/T_Tc_yreal(i),2.0)
    delta_modification2(i)=delta_T_Tc(i).*T_Tc_xreal_inv(i)
    delta_modification(i)=(delta_T_Tc(i)).*sqrt(T_Tc_xreal_inv(i))
end                                  
                                     
assign_a= -0.63             
assign_b= 0.99             
assign_c= 0.97             
assign_d= -0.63             
assign_e= 0.28             
                                     
figure(1)                        
i=1:length(T_Tc_xreal);              
T_Tc_London_Modified(i)=1.0./sqrt(1.0-assign_a*T_Tc_xreal5(i)-(assign_b)*T_Tc_xreal4(i)-(assign_c)*T_Tc_xreal3(i)-(assign_d)*T_Tc_xreal2(i)-(assign_e)*T_Tc_xreal(i));
plot(T_Tc_xreal(i), T_Tc_yreal(i),'ko');
hold on
plot(T_Tc_xreal(i), T_Tc_London_Modified(i), 'k-');
hold on;
xlabel('T/Tc');
ylabel('Lambda(T)/Lambda(0)');
axis([0.20 1.00 0.90 3.50])

figure(2)                        
plot(T_Tc_xreal(i), T_Tc_yreal(i),'ko');
hold on
plot(T_Tc_xreal(i), T_Tc_London_Modified(i), 'k-');
hold on;
xlabel('T/Tc');
ylabel('Lambda(T)/Lambda(0)');
axis([0.20 1.00 0.90 3.50])

figure(3)                        
i=1:length(T_Tc_xreal);
plot(T_Tc_xreal(i), T_Tc_yreal_inv(i), 'ko');
hold on;
plot(T_Tc_xreal(i), 1./T_Tc_London_Modified(i), 'k-');
hold on;
xlabel('T/Tc');
ylabel('Lambda(0)/Lambda(T)');
axis([0.20 1.00 0.20 1.01])

figure(4)                        
i1=1:1:100;                          
i002(i1)=0.01*i1;                    
me_m0(i1)=power(1./(i002(i1).*i002(i1))-(assign_a)*power(i002(i1),3.0)-(assign_b).*power(i002(i1),2.0)-(assign_c)*power(i002(i1),1.0)-(assign_d)-(assign_e)*power(i002(i1),-1.0), 2.0);
i2=1:1:100;                          
i002(i2+100)=0.01*i2+1.00;           
me_m0(i2+100)=power(1.0-(assign_a)*power(i002(i2+80),-5.0)-(assign_b)*power(i002(i2+80),-4.0)-(assign_c)*power(i002(i2+80),-3.0)-(assign_d)*power(i002(i2+80),-2.0)-(assign_e)*power(i002(i2+80),-1.0), 2.0);
i=1:1:200;                           
plot(i002(i), me_m0(i), 'k-');       
hold on;                             
xlabel('T/Tc');                      
ylabel('Me/M0');                     
axis([0 2.00 -10.0 1000])             

figure(5)                        
plot(i002(i), log10(me_m0(i)), 'k-');
hold on;
xlabel('T/Tc');
ylabel('Log10(Me/M0)')
axis([0 2.00 -4.00 4.00])
