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);