1clear;
2clc;
3
4initialDeveloperConcentration = single(1);
5reservoirThickness = single(1000);
6activeLayerThickness = single(0.1);
7crystalsPerPixel = single(500);
8initialCrystalRadius = single(0.00001);
9initialSilverSaltDensity = single(1);
10developerConsumptionConst = single(2000000);
11crystalGrowthConst = single(0.00001);
12silverSaltConsumptionConst = single(2000000);
13totalDevelopmentTime = single(100);
14agitateCount = single(1);
15developmentSteps = single(12);
16filmArea = single(864);
17sigmaConst = single(0.2);
18layerMixConst = single(0.5);
19layerTimeDivisor = single(20);
20
21numAmps = 5;
22numMixes = 10;
23
24sameChannelPeak = zeros(numAmps,numMixes);
25otherChannelPeak = zeros(numAmps,numMixes);
26curveProfile = zeros(1000,numAmps,numMixes);
27
28pulseAmplitudes = (2*ones(1,numAmps)).^linspace(6,17,numAmps);
29layerMixConsts = linspace(0,1,numMixes).^2;
30layerMixCoefs = layerMixConsts.^(totalDevelopmentTime/layerTimeDivisor);
31for mix = 1:numMixes
32    for amp = 1:numAmps
33      thisAmplitude = pulseAmplitudes(amp);
34      thisLayerMixConst = layerMixConsts(mix);
35
36      initialData = single(zeros(1000,1000,10));
37      initialData(:,:,1:3) = initialCrystalRadius;
38      initialData(:,:,4:6) = 0*crystalsPerPixel*0.00015387105;
39      initialData(500,:,4) = thisAmplitude*crystalsPerPixel*0.00015387105;
40      initialData(:,:,7:9) = initialSilverSaltDensity;
41      initialData(:,:,10)  = initialDeveloperConcentration;
42      reservoirConcentration = initialDeveloperConcentration;
43
44      for i = 1:developmentSteps
45         outData = single(zeros(1000,1000,10));
46         outReservoirConcentration = single(ones(2,1));
47         filmulateIterationGenerator(reservoirConcentration,reservoirThickness, ...
48                                     crystalGrowthConst,activeLayerThickness, ...
49                                     developerConsumptionConst,silverSaltConsumptionConst, ...
50                                     totalDevelopmentTime/developmentSteps,filmArea,sigmaConst, ...
51                                     thisLayerMixConst,layerTimeDivisor,true, ...
52                                     initialData,outData,outReservoirConcentration);
53         initialData = outData;
54         reservoirConcentration = outReservoirConcentration(1);
55      end
56
57      sameChannelPeak(amp,mix) = outData(500,500,1);
58      otherChannelPeak(amp,mix) = mean(outData(500,500,2:3));
59      curveProfile(:,amp,mix) = outData(:,500,2);
60    end
61end
62
63
64% figure(1);
65% imagesc(outData(:,:,2));
66%
67% xPts = -499:500;
68%
69% figure(2);
70% baseline = outData(1,500,2);
71% peak = outData(500,500,2);
72% sigma = 30;
73% gaussY = baseline+(peak-baseline)*exp(-xPts.^2/(2*sigma.^2));
74% plot(xPts,outData(:,500,2),'b',xPts,gaussY,'r')
75
76figure(2);
77mesh(squeeze(curveProfile(1000,:,:)))
78title('0 Points');
79
80figure(3);
81% plot(pulseAmplitudes,sameChannelPeak-otherChannelPeak,'b',pulseAmplitudes,0.001./(1+(pulseAmplitudes/100000).^(0.75))-0.001,'r');
82% plot(layerMixConsts,sameChannelPeak-otherChannelPeak,'b',layerMixConsts, 0.001./(1 + (thisAmplitude/100000).^(0.75) + (layerMixConsts/5000).^(0.3))-0.001,'r');
83mesh(sameChannelPeak - squeeze(curveProfile(1000,:,:)))
84title('Same channel additional peak');
85
86baseline = 0.00101;
87%baseline = squeeze(curveProfile(1000,:,:))
88figure(4);
89% plot(pulseAmplitudes,otherChannelPeak,'b',pulseAmplitudes,0.00101./(1+(pulseAmplitudes/2000000000).^(0.4)),'r');
90mesh((otherChannelPeak - baseline)/(sameChannelPeak(:,2:end) - baseline))
91title('Other channel peak ratio');
92
93figure(5);
94% plot(pulseAmplitudes,otherChannelPeak,'b',pulseAmplitudes,0.00101./(1+(pulseAmplitudes/2000000000).^(0.4)),'r');
95mesh(layerMixCoefs, log2(pulseAmplitudes), otherChannelPeak)
96title('Other channel peak ratio');
97
98figure(6);
99plot(squeeze(curveProfile(:,end-1,:)));
100
101figure(7);
102plot(layerMixCoefs,otherChannelPeak(end,:))
103
104% figure(6);
105% plot(sameChannelPeak,otherChannelPeak);