%Phoebe Fogelman
% 27 March 2015
% University of Tennessee, Knoxville - Sharpe Lab
% Move Excel calculations to matlab
%frames=[60,65,70,75,78,80,83,85,88,90,92,95,98,100];
%frames=[9 10 12 14 16 20 22 24 26 28 32 34 38];
%%
%Get user input
unique_name=inputdlg('Enter job name or alternate analysis identifier.');
unique_name=unique_name{1};
%Define Material Properties 
mat_name='A508_neg85C';
E=212100; %N/mm2
nu=0.3;
s_y=555;%N/mm^2
n=8;%exponent
s1=-.11111;
s2=0.07089;
s3=0.25289;	
sigma1=2.418;
sigma2=0.3186;
sigma3=-5.6213;
In=4.67567;
L=1;
j_mat=load(sprintf('%s_j_integral.dat',unique_name));
%j_mat=dlmread('j_integral.txt','',2,0);
r_1mm=j_mat/s_y;
r_5mm=5*r_1mm;

%constraints for getting A2's to average
r_min=1;
r_max=4;

frames=dlmread(sprintf('%s_MATLAB_PARAMS.txt',unique_name),',',1,0);
frames=frames(1:end-1);
num_frames=length(frames);

a_J_mat=zeros([num_frames,2]);
use_pts=0;
for i=1:num_frames
     fname=sprintf('Frame%d',frames(i));
    stress_file=sprintf('%s_frame%drpt.txt',unique_name,frames(i));
    stress_table=readtable(stress_file,'Headerlines',3,'ReadVariableNames',0);
    stress_cell=table2cell(stress_table);
    x_pts=length(stress_cell);
    calcs.(fname).stress_path=zeros([x_pts, 2]);
    for q=1:x_pts
        stress_cell{q}=strsplit(stress_cell{q},' ');
        inter_mat=str2double(stress_cell{q});
        calcs.(fname).stress_path(q,1)=inter_mat(1);
        calcs.(fname).stress_path(q,2)=inter_mat(2);
    end
    %r*yield strength/J
    calcs.(fname).R_J=calcs.(fname).stress_path(:,1)./j_mat(frames(i),2).*s_y;
    %point stress/yield 
    x_pts=length(calcs.(fname).R_J);
    calcs.(fname).ptstress_Ys=calcs.(fname).stress_path(:,2)./s_y;
    %A2^2 coefficient = (radius/L)^s3*sigma3
    calcs.(fname).A2_2_coe=(calcs.(fname).stress_path(:,1)./L).^s3.*sigma3;
    %other two coefficients using s2,sigma2 and s3,sigma3 in place of
    %s1,sigma1
    calcs.(fname).A2_1_coe=(calcs.(fname).stress_path(:,1)./L).^s2.*sigma2;
    calcs.(fname).A2_0_coe=(calcs.(fname).stress_path(:,1)./L).^s1.*sigma1;
    calcs.(fname).A1=(j_mat(frames(i),2)./((s_y/E)*s_y*In*L))^(-1*s1);
    %sigma/y_s constant
    calcs.(fname).sig_Y_s=calcs.(fname).ptstress_Ys./calcs.(fname).A1;
    %quadratic_addition
    calcs.(fname).quad_add=(-1*calcs.(fname).A2_1_coe+sqrt( calcs.(fname).A2_1_coe.^2-(4* calcs.(fname).A2_2_coe.*( calcs.(fname).A2_0_coe- calcs.(fname).sig_Y_s))))./(2* calcs.(fname).A2_2_coe);
    %quadratic_subtraction
    calcs.(fname).quad_subtr=(-1* calcs.(fname).A2_1_coe-sqrt( calcs.(fname).A2_1_coe.^2-(4* calcs.(fname).A2_2_coe.*( calcs.(fname).A2_0_coe- calcs.(fname).sig_Y_s))))./(2* calcs.(fname).A2_2_coe);
    calcs.(fname).pass_A2s=zeros([x_pts 1]);
    count=0;
    for j=1:x_pts
        if calcs.(fname).R_J(j)>=2 && calcs.(fname).R_J(j)<=5
            count=count+1;
            calcs.(fname).pass_A2s(count)=calcs.(fname).quad_add(j);
        end
    end
    calcs.(fname).pass_A2s=calcs.(fname).pass_A2s(1:count);
    calcs.(fname).A2_avg=mean(abs(calcs.(fname).pass_A2s));
    if calcs.(fname).A2_avg<=1 && j_mat(frames(i),2)<=100 && calcs.(fname).A2_avg>=0.05
        use_pts=use_pts+1;
        a_J_mat(use_pts,1)=j_mat(frames(i),2); 
        a_J_mat(use_pts,2)=calcs.(fname).A2_avg ;
    end
end
a_J_mat=a_J_mat(1:use_pts,:);
plot(a_J_mat(:,2),a_J_mat(:,1));
xlabel('|A2|');
ylabel('J (N/mm)');
xlim([0 .45]);ylim([0 140]);

%%
%Experimental Material Failure Curve
%Define critical stress value - determined from Chao RPV Paper based on
%results from Hohe paper for 3 point bending 
mat_fx=[0.15
0.151
0.152
0.153
0.154
0.155
0.156
0.157
0.158
0.159
0.16
0.161
0.162
0.163
0.164
0.165
0.166
0.167
0.168
0.169
0.17
0.171
0.172
0.173
0.174
0.175
0.176
0.177
0.178
0.179
0.18
0.181
0.182
0.183
0.184
0.185
0.186
0.187
0.188
0.189
0.19
0.191
0.192
0.193
0.194
0.195
0.196
0.197
0.198
0.199
0.2
0.201
0.202
0.203
0.204
0.205
0.206
0.207
0.208
0.209
0.21
0.211
0.212
0.213
0.214
0.215
0.216
0.217
0.218
0.219
0.22
0.221
0.222
0.223
0.224
0.225
0.226
0.227
0.228
0.229
0.23
0.231
0.232
0.233
0.234
0.235
0.236
0.237
0.238
0.239
0.24
0.241
0.242
0.243
0.244
0.245
0.246
0.247
0.248
0.249
0.25
0.251
0.252
0.253
0.254
0.255
0.256
0.257
0.258
0.259
0.26
0.261
0.262
0.263
0.264
0.265
0.266
0.267
0.268
0.269
0.27
0.271
0.272
0.273
0.274
0.275
0.276
0.277
0.278
0.279
0.28
0.281
0.282
0.283
0.284
0.285
0.286
0.287
0.288
0.289
0.29
0.291
0.292
0.293
0.294
0.295
0.296
0.297
0.298
0.299
0.3
0.301
0.302
0.303
0.304
0.305
0.306
0.307
0.308
0.309
0.31
0.311
0.312
0.313
0.314
0.315
0.316
0.317
0.318
0.319
0.32
0.321
0.322
0.323
0.324
0.325
0.326
0.327
0.328
0.329
0.33
0.331
0.332
0.333
0.334
0.335
0.336
0.337
0.338
0.339
0.34
0.341
0.342
0.343
0.344
0.345
0.346
0.347
0.348
0.349
0.35
0.351
0.352
0.353
0.354
0.355
0.356
0.357
0.358
0.359
0.36
0.361
0.362
0.363
0.364
0.365
0.366
0.367
0.368
0.369
0.37
0.371
0.372
0.373
0.374
0.375
0.376
0.377
0.378
0.379
0.38
0.381
0.382
0.383
0.384
0.385
0.386
0.387
0.388
0.389
0.39
0.391
0.392
0.393
0.394
0.395
0.396
0.397
0.398
0.399
0.4
0.401
0.402
0.403
0.404
0.405
0.406
0.407
0.408
0.409
0.41
0.411
0.412
0.413
0.414
0.415
0.416
0.417
0.418
0.419
0.42
0.421
0.422
0.423
0.424
0.425
0.426
0.427
0.428
0.429
0.43
0.431
0.432
0.433
0.434
0.435
0.436
0.437
0.438
0.439
0.44
0.441
0.442
0.443
0.444
0.445
0.446
0.447
0.448
0.449
0.45
0.451
0.452
0.453
0.454
0.455
0.456
0.457
0.458
0.459];
mat_fy=[11.31790647
11.35977982
11.40204343
11.44470079
11.48775543
11.53121092
11.57507086
11.61933893
11.66401883
11.7091143
11.75462913
11.80056718
11.84693232
11.8937285
11.94095969
11.98862993
12.03674331
12.08530396
12.13431605
12.18378383
12.23371159
12.28410367
12.33496445
12.3862984
12.43811001
12.49040385
12.54318453
12.59645672
12.65022517
12.70449465
12.75927003
12.81455621
12.87035817
12.92668094
12.98352962
13.04090938
13.09882544
13.15728309
13.21628769
13.27584468
13.33595954
13.39663784
13.45788522
13.51970739
13.58211011
13.64509926
13.70868074
13.77286057
13.83764482
13.90303965
13.96905129
14.03568607
14.10295036
14.17085066
14.23939352
14.30858559
14.37843359
14.44894435
14.52012477
14.59198184
14.66452265
14.73775438
14.81168429
14.88631976
14.96166824
15.03773728
15.11453455
15.19206781
15.2703449
15.34937379
15.42916256
15.50971936
15.59105247
15.6731703
15.75608132
15.83979417
15.92431756
16.00966033
16.09583145
16.18283999
16.27069515
16.35940626
16.44898275
16.53943421
16.63077034
16.72300097
16.81613605
16.9101857
17.00516015
17.10106976
17.19792506
17.2957367
17.39451549
17.49427236
17.59501842
17.69676492
17.79952326
17.90330499
18.00812185
18.11398569
18.22090858
18.32890271
18.43798046
18.54815439
18.6594372
18.77184181
18.88538129
19.00006889
19.11591807
19.23294245
19.35115587
19.47057233
19.59120605
19.71307145
19.83618314
19.96055595
20.08620492
20.21314528
20.3413925
20.47096227
20.60187048
20.73413328
20.86776703
21.00278832
21.139214
21.27706113
21.41634705
21.55708933
21.69930579
21.84301452
21.98823387
22.13498245
22.28327915
22.43314313
22.58459383
22.73765097
22.89233456
23.04866493
23.20666266
23.36634867
23.52774419
23.69087073
23.85575016
24.02240465
24.19085671
24.36112917
24.53324521
24.70722838
24.88310255
25.06089195
25.2406212
25.42231526
25.60599949
25.79169963
25.9794418
26.16925253
26.36115873
26.55518775
26.75136734
26.94972567
27.15029136
27.35309345
27.55816145
27.7655253
27.97521541
28.18726267
28.40169844
28.61855457
28.83786341
29.05965781
29.28397112
29.51083724
29.74029058
29.9723661
30.20709931
30.44452628
30.68468364
30.92760862
31.17333902
31.42191326
31.67337035
31.92774995
32.18509233
32.44543842
32.7088298
32.97530873
33.24491813
33.51770164
33.79370359
34.07296902
34.35554374
34.64147425
34.93080786
35.22359263
35.51987741
35.81971184
36.12314641
36.4302324
36.74102197
37.05556813
37.37392477
37.69614668
38.02228955
38.35241002
38.68656566
39.024815
39.36721758
39.71383391
40.06472554
40.41995506
40.7795861
41.1436834
41.51231277
41.88554117
42.26343667
42.64606854
43.03350721
43.42582433
43.82309278
44.22538671
44.63278154
45.045354
45.46318214
45.8863454
46.31492456
46.74900184
47.18866091
47.63398688
48.08506637
48.54198752
49.00484004
49.47371521
49.94870594
50.4299068
50.91741401
51.41132555
51.91174112
52.41876223
52.93249219
53.45303619
53.98050131
54.51499656
55.05663294
55.60552345
56.16178316
56.72552923
57.29688097
57.87595986
58.46288962
59.05779625
59.66080808
60.27205579
60.89167249
61.51979377
62.15655773
62.80210506
63.45657906
64.12012572
64.79289379
65.47503479
66.16670311
66.86805605
67.57925391
68.30046002
69.03184079
69.77356586
70.52580806
71.28874356
72.0625519
72.84741608
73.64352262
74.45106166
75.27022702
76.10121627
76.94423084
77.79947608
78.66716136
79.54750015
80.44071012
81.34701322
82.26663577
83.19980858
84.14676703
85.10775118
86.08300588
87.07278084
88.0773308
89.09691558
90.13180025
91.18225519
92.24855628
93.33098495
94.42982836
95.54537951
96.67793737
97.82780702
98.99529979
100.1807334
101.3844322
102.606727
103.8479556
105.108463
106.388601
107.688729
109.0092137
110.3504296
111.7127592
113.0965927
114.5023288
115.9303745
];
hold on 
plot (mat_fx(:),mat_fy(:),'r')
%%
%Curve fit Data
cf=questdlg('Curve fit and add experimental data to graph?');
if strcmp(cf,'Yes')
close all
plot(a_J_mat(:,2),a_J_mat(:,1),'k*','MarkerSize',3);
hold on
plot (mat_fx(:),mat_fy(:),'r');
best_fit=polyfit(a_J_mat(:,2),a_J_mat(:,1),4);
x=0.29:0.01:0.36; y=polyval(best_fit,x);plot(x,y,'b');
%Add experimental Data -- %this is for deep 3PB - CHANGE to reflect a
%different specimen
exp_a2=0.292.*ones([5,1]);
exp_j=[45.96014851
43.42068647
51.91419142
31.87977558
82.4189703];
plot(exp_a2,exp_j,'o','MarkerFaceColor',[.87 0 .35], 'MarkerEdgeColor',[.87,0,.35],'MarkerSize',3);
xlabel('|A2|');
ylabel('J (N/mm)');
xlim([0 .45]);ylim([0 140]);
legend('Abaqus Output','Material Failure Curve','Best Fit Line','Experimental Data');
title('J-A2 Plot for Deep-Cracked A508 Specimen at -85 deg. C')
end