Problem with script in if -else conditions

I am PhD student. I have to calculate a parameter (settling velocity-ws ) using two different sets of equations selected according a criteria (sediment diameter - D50). Do not worry about the subject. But the idea is that ws depends on D50. For sand (D50>0.00005), I want to use one set of equations to calculate ws. For other sediment particles (D50<0.00005), I want to use another set of equations to calculate ws. I load all sediment diameter (D50) values using a text file.

I used if and else conditions (please see the script written by me; I am new to matlab, this may not be efficient).

However, the results from my script gives an answer but using only one set of equations to calculate settling velocity. There should be a problem of writing if condition in my script. Here is the sample data set.

D50=
0.0002626
0.0002626
0.000003504
0.0000108
0.0000985

This is the answer;

ws=
0.069
0.069
0.0000123
0.000117
0.009717

The correct answer should be;

ws=
0.03681
0.03681
0.0000123
0.000117
0.009717

This is my code for actual set of equations; I want the answers written into a text file.

    %%%Settling velocity

g=9.81;

psed=2650;

pw=1027;

neu=0.000001004;

p=0.4;

load D50.txt

          if D50 > 0.0000625
          Dstar=D50*(g*((psed/pw)-1)/neu^2)^(1/3);
          Dstarpower=Dstar.^3;  
          idleA=ones(size(D50));   
          idleB=107.3296*idleA+1.049*Dstarpower; 
          idleC=idleB.^0.5;
          idleD=idleC-10.36*idleA; 
          idleE=neu*idleD;
          ws=idleE./D50

      else 
          ws=1000000*(D50.^2)
      end

      CoefA=ws/((1-p)*psed);

      fileID1=fopen('ws.txt','w');
      fprintf(fileID1,'%6.6f\r\n',ws);
      fclose(fileID1);

fileID2=fopen('CoefA.txt','w');
fprintf(fileID2,'%6.6f\r\n',CoefA);
fclose(fileID2);

I am very grateful if you can sort out this problem. Do not bother about equations, if equations are complex please give me a sample code for a case

if D50>0.00005
    D_STAR=100*D50;
    ws=(D_STAR*100+10)/.D50
else
     ws=1000000*(D50.^2)
end

Thank you. Ruwan Sampath

Answer 1

[Andrew Newell]

The problem is that you are testing the entire array at once with your code. If you type

D50 > 0.0000625

you get

ans =

     1
     1
     0
     0
     1

An if test is only true if all the elements are true. You need to create a loop to test each element separately:

ws = zeros(size(D50));
for ii = 1:numel(D50)
  if false %D50(ii) > 0.0000625
      Dstar=D50(ii)*(g*((psed/pw)-1)/neu^2)^(1/3);
      Dstarpower=Dstar.^3;
      idleA=ones(size(D50(ii)));
      idleB=107.3296*idleA+1.049*Dstarpower;
      idleC=idleB.^0.5;
      idleD=idleC-10.36*idleA;
      idleE=neu*idleD;
      ws(ii)=idleE./D50(ii);
  else
      ws(ii)=1000000*(D50(ii).^2);
  end
end
ws

The result is

ws =

   0.036807126831593
   0.036807126831593
   0.000012278016000
   0.000116640000000
   0.007330330176007

The last entry disagrees with your expectation for the correct answer, but it is consistent with the test condition.

[Sampath Dissanayake]

Thank you very much for your answer. It worked well.

Answer 2

[Star Strider]

I believe some of your ‘correct’ answers aren’t actually correct. The else block applies to some of the ‘D50’ elements.

The easiest solution (without disrupting your code because I don’t understand what you’re doing) is to use a for loop to iterate through your ‘D50’ vector:

for k1 = 1:length(D50)
    D50v = D50(k1);
          if D50v > 0.0000625
          Dstar=D50v*(g*((psed/pw)-1)/neu^2)^(1/3);
          Dstarpower=Dstar.^3;  
          idleA=ones(size(D50v));   
          idleB=107.3296*idleA+1.049*Dstarpower; 
          idleC=idleB.^0.5;
          idleD=idleC-10.36*idleA; 
          idleE=neu*idleD;
          ws(k1)=idleE./D50v;

        else 
            ws(k1)=1000000*(D50v.^2);
        end
  end

This produces for ‘ws’:

ws =

      36.8071e-003    36.8071e-003    12.2780e-006   116.6400e-006     7.3303e-003

[Sampath Dissanayake]

Thank you very much for your answer. It worked well.

[Star Strider]

My pleasure.