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

 clear all;                          
                                     
T_Tc_inv_x00=149;                    
T_Tc_inv_x02=254;                    
T_Tc_inv_x04=359;                    
T_Tc_inv_x06=465;                    
T_Tc_inv_x08=570;                    
T_Tc_inv_x10=675;                    
T_Tc_inv_y00=407;                    
T_Tc_inv_y02=331;                    
T_Tc_inv_y04=256;                    
T_Tc_inv_y06=181;                    
T_Tc_inv_y08=105;                    
T_Tc_inv_y10=30;                     
                                     
T_Tc_inv_x=[175 182 189 196 203 209 216 223 230 236 243 250 257 263 270 277 284 290 297 304 311 318 324 331 337 345 351 358 365 372 378 385 392 398 405 412 419 425 432 439 446 452 459 466 473 480 486 493 500 506 513 520 527 533 540 547 553 561 567 574 581 588 594 601 608 615 621 628 635 641 648 655 662 668 675];
T_Tc_inv_y=[32 33 35 36 39 41 43 46 49 52 55 59 63 67 71 75 81 86 91 95 101 107 112 118 125 130 137 144 150 156 163 170 178 184 191 199 206 214 221 229 236 243 250 258 264 271 279 285 291 299 304 310 318 324 329 335 339 347 350 355 360 365 370 372 376 380 384 387 391 394 396 399 401 403 405];
                                     
x_inv_01=(T_Tc_inv_x10-T_Tc_inv_x00)/10.0
y_inv_01=(T_Tc_inv_y00-T_Tc_inv_y10)/10.0
                                     
vector_size=0;                       
                                     
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_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.341             
assign_b= -0.931             
assign_c= 0.773             
assign_d= 1.279             
assign_e= 0.183             
                                     
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), 1.0./power(T_Tc_yreal(i),2.0),'ko');
hold on
plot(T_Tc_xreal(i), 1.0./power(T_Tc_London_Modified(i),2.0), 'k-');
hold on;
xlabel('T/Tc');
ylabel('Lambda(T)/Lambda(0)');
axis([0 1.00 0.0 1.00])

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 1.00 0.90 10.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 1.00 0 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])
