invreyfunc.m

function answer=invreyfunc(y,Rey,H)

if Rey~=-1      %A CONSTANT TURBULENT REYNOLDS NUMBER CAN BE SIMULATED. SIMPLY PUT Rey~=-1 IN linearstability.m
    answer(1,1)=1/Rey;
    answer(1,2)=0;
    answer(1,3)=0;
elseif Rey==-1
    
    h=0.69;
    
    %TURBULENT VISCOSITY PROFILE FROM THE LES
    %HEIGH MADE DIMENSIONLESS USING CANOPY HEIGHT
    datay=1/h*[0.0400200001900000,0.1199909970000000,0.2000309974000000,0.2800019979000000,0.3599730134000000,0.4400130212000000,0.5199840069000000,0.6000239849000000,0.6799950004000000,0.7599660158000000,0.8400059938000000,0.9199770093000000,1.0000170470000000,1.0799880030000000,1.1600279810000000,1.2399989370000000,1.3199700120000000,1.4000099900000000,1.4799810650000000,1.5600210430000000,1.6399919990000000,1.7200319770000000,1.8000030520000000,1.8799740080000000,1.9600139860000000,2.0399849410000000,2.1200249200000000,2.1999959950000000,2.2799670700000000,2.3600070480000000,2.4399781230000000,2.5200181010000000,2.5999889370000000,2.6800289150000000,2.7599999900000000,2.8399710660000000,2.9200110440000000,2.9999818800000000,3.0800218580000000,3.1599929330000000,3.2400331500000000,3.3200042250000000,3.3999748230000000,3.4993348120000000,3.6333329680000000,3.7943100930000000,3.9756419660000000,4.1716022490000000,4.3772912030000000,4.5890522000000000,4.8039178850000000,5.0200948720000000];
    %KINEMATIC TURBULENT VISCOSITY PROFILE OBTAINED AT Ucanopy_top=3.8 m/s.
    %MADE INTO THE INVERSE OF A REYNOLDS NUMBER USING Ucanopy_top AND h.
    %THE TURBULENT REYNOLDS NUMBER IS CONSTANT WITH FLOW VELOCITY IN THE
    %LES, SO WE TAKE IT CONSTANT HERE IN THE LINEAR ANALYSIS.
    datainvrey=1/h/3.8*[0.0008077020175000,0.0016784397890000,0.0020962969870000,0.0025884245990000,0.0031854400880000,0.0039366465060000,0.0049014044930000,0.0059775398110000,0.0072549227630000,0.0082411756740000,0.0081031452860000,0.0081592313950000,0.0080548273400000,0.0080038765450000,0.0079035898670000,0.0078218150880000,0.0077148280110000,0.0076171872210000,0.0075042941610000,0.0073980130260000,0.0072838715280000,0.0071822293100000,0.0070902747100000,0.0069853612220000,0.0068887188100000,0.0067939707080000,0.0067016836260000,0.0066129318440000,0.0065276199020000,0.0064415135420000,0.0063616726550000,0.0062879375180000,0.0062128449790000,0.0061428295450000,0.0060881599780000,0.0060402243400000,0.0059983148240000,0.0059896460730000,0.0059726093900000,0.0060217198920000,0.0060646045020000,0.0062577901410000,0.0064180730840000,0.0072373519650000,0.0079143466430000,0.0084540322420000,0.0088639426980000,0.0091571537780000,0.0093660280110000,0.0095055401330000,0.0096066938710000,0.0096680037680000];

    answer(1,1)=pchip(datay,datainvrey,y);      %INVERSE OF THE REYNOLDS NUMBER

    hdata=max(datay);                           %MAXIMUM VALUE IN y WHERE THE REYNOLDS PROFILE IS DEFINED

    %THE VALUE IS INTERPOLATED USING THE pchip METHOD OF MATLAB. OUTSIDE
    %THE RANGE OF y WHERE THE PROFILE IS DEFINED, WE USE A SIMPLE LINEAR
    %EXTRAPOLATION
    if y<hdata
        imethod='pchip';
    else
        imethod='linear';
    end

    answer(1,1)=interp1(datay,datainvrey,y,imethod,'extrap');
    delta=(0.1199909970000000-0.0400200001900000)/h;
    uc=answer(1,1);
    up=interp1(datay,datainvrey,y+delta,imethod,'extrap');
    if y<delta
        um=0;
    else
        um=interp1(datay,datainvrey,y-delta,imethod,'extrap');
    end
    answer(1,2)=(up-um)/2/delta;                %FIRST DERIVATIVE OF THE INVERSE OF THE REYNOLDS NUMBER
    answer(1,3)=(up-2*uc+um)/delta^2;           %SECOND DERIVATIVE OF THE INVERSE OF THE REYNOLDS NUMBER
else
    bigerror=invalidvalueofReynolds
    return
end
return