function AREA_TEMPLATE = rbcA0(ipts) % Calulates the dimensioless area of a resting red cell % as a function of radial position % The range of r0 is zero to one. % The integration is divided up into five regions with different spacings. % Each section is divided into ipts, % Near the origin and the rim, the spacings of r0 values are small because of singularities. % The last interval stops just short of r0 = 1, and % the last point at r0 = 1 is calculated by multiplying % the height of the last point times the radius (1.0) and adding % that area to the preceding value. %AREA_TEMPLATE is a four row by (5*ipts+1) matrix. The first row contains % the dimensionless r0, the second row contains the dimensionless area, the % third row contains "s" and the last row contains theta. % The total interval is divided into five regions each have "ipts" elements. r0 = linspace(0,.1-.1/ipts, ipts); % Each interval stops just short of the starting r1 = linspace(0.1,.9-.9/ipts,ipts); % value for the next interval. r2 = linspace(.9,.99-.09/ipts,ipts); r3 = linspace(.99,.999-.009/ipts,ipts); r4 = linspace(.999,.9999,ipts); % First zero out the global array so we don't have leftover values from a previous try. AREA_TEMPLATE = [0; 0; 0; 0; 0]; for i = 1:ipts if (i==1) % the first point in the first region is zero, so just set A0 = 0. A0 =0; S0 = 0; Theta0 = 0; z0 = rbcth(0); else A0 = integral(@dA0dr0,0,r0(i)); %Otherwise use the Matlab function quadl to calculate S0 = integral(@dsdr,0,r0(i)); %the area and s from zero up to the current (ith) value of r0 Theta0 = atan(-rbcdzdr(r0(i))); z0 = rbcth(r0(i)); end A1 = integral(@dA0dr0,0,r1(i)); % for each of the five intervals S1 = integral(@dsdr,0,r1(i)); Theta1 = atan(-rbcdzdr(r1(i))); z1 = rbcth(r1(i)); A2 = integral(@dA0dr0,0,r2(i)); S2 = integral(@dsdr,0,r2(i)); Theta2 = atan(-rbcdzdr(r2(i))); z2 = rbcth(r2(i)); A3 = integral(@dA0dr0,0,r3(i)); S3 = integral(@dsdr,0,r3(i)); Theta3 = atan(-rbcdzdr(r3(i))); z3 = rbcth(r3(i)); A4 = integral(@dA0dr0,0,r4(i)); S4 = integral(@dsdr,0,r4(i)); Theta4 = atan(-rbcdzdr(r4(i))); z4 = rbcth(r4(i)); % Fill the area Template with r0 in the first row and A0 in the second for the ith value of each region AREA_TEMPLATE(:,i) = [r0(i); A0; S0; Theta0; z0]; % Assign the values for the first region (0 to just short of .1)) AREA_TEMPLATE(:,i+ipts) = [r1(i); A1; S1; Theta1; z1]; % The second region (.1 just short of .9) starts at column (ipts + 1) AREA_TEMPLATE(:,(i+2*ipts)) = [r2(i); A2; S2; Theta2; z2]; % The third region (.9 to just short of .99) starts at 2*ipts + 1 AREA_TEMPLATE(:,(i+3*ipts)) = [r3(i); A3; S3; Theta3; z3]; % Fourth region (.99 to just short of .999) starts at 3*ipts + 1 AREA_TEMPLATE(:,(i+4*ipts)) = [r4(i); A4; S4; Theta4; z4]; % Fifth region (.999 to just short of .9999) starts at 4*ipts + 1 end % fill the last array column by a trapezoid step AREA_TEMPLATE(1,5*ipts+1) = 1; % The last value of r0 is one. % The last area increment is the thickness at the last point in region 5 times the radius, % which is the value half way between the last value of the fifth interval, and the end value (1.0). AREA_TEMPLATE(2,5*ipts+1) = AREA_TEMPLATE(2,5*ipts) + rbcth(r4(ipts))*(1+r4(ipts))/2; AREA_TEMPLATE(3,5*ipts+1) = AREA_TEMPLATE(3,5*ipts) + rbcth(r4(ipts)); AREA_TEMPLATE(4,5*ipts+1) = pi/2; %Last value of theta is pi/2 AREA_TEMPLATE(5,5*ipts+1) = 0; %Last value of z is zero disp('Area Template is set.'); format shortG; % disp( AREA_TEMPLATE');