%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%code name: feasibility.cc

%manuel (Manuel.Lonfat@rms.com)
%last modified: 01.14.2005

%Test code for feasibility study of rainfall modeling. R-CLIPER-like, plus 
%contributions from Fourier orders 1 and 2.

%Inputs:
%. storm track, intensity, size and speed info at 6-hour resolution, with info on
%  amplitude and direction of the environmental shear.
%. file with info on mean distribution of rain in TCs of various intensity.
%. files with Fourier coefficients, first and second order, shear-relative.

%Outputs:
%. array of accumulated rainfall on a 1/20-degree resolution, over the entire 
%  storm history.
%. same but at hourly resolution: hourly increments in rainfall inputs.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Construction of the master grid, at 1/10-degree resolution,   %
% defined between 50,100 West, and 10, 60 N                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

for i=1:50*10,
  for j=1:50*10,
    masterLat(i,j) = 10.0+j/10.0;
    masterLon(i,j) = -100.0+i/10.0;
    master_history(i,j)=0.0;
  end
end

% Track data for Hurricane Charley, 2004                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

track_data=load('frances2004.txt');
[mm,nn]=size(track_data);

% data interval (in hours)
interv       = 6;
% integration time step (in fraction of hours)
timeStep     = 0.25;
% number of intervals between 2 data points
% 48 if timeStep=0.25 and interv=12
numbOfInterv = (1/timeStep)*interv;


% Loop over times within the track file                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  

%range of interest (in input file) for land accumulation, roughly...
%charley - init=10 to mm=eof/25
%frances - init=33 to mm=eof/68
%ivan    - init=52 to mm=eof/87 (64 for first landmass)
%jeanne  - init=35 to mm=eof/64
%allison - init=1 to mm=eof/25

init=1;
mm=68;
index=1;
for m=init:mm-1,
  m
  interpTime      = 0;
  interpLat       = track_data(m,6);
  interpLon       = track_data(m,7);
  interpShearAmpl = track_data(m,10);
  interpShearDir  = track_data(m,11);
  interpRmax      = track_data(m,12);
  
  while interpTime<interv,
  
% Linear interpolation of track database: 15-minute resolution         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    interpLat       = track_data(m,6)+(track_data(m+1,6)-track_data(m,6))*((1/timeStep)*interpTime)/numbOfInterv;
    interpLon       = track_data(m,7)+(track_data(m+1,7)-track_data(m,7))*((1/timeStep)*interpTime)/numbOfInterv;
    interpShearAmpl = track_data(m,10)+(track_data(m+1,10)-track_data(m,10))*((1/timeStep)*interpTime)/numbOfInterv;
    interpWind      = track_data(m,9)+(track_data(m+1,9)-track_data(m,9))*((1/timeStep)*interpTime)/numbOfInterv; 
    interpRmax      = track_data(m,12)+(track_data(m+1,12)-track_data(m,12))*((1/timeStep)*interpTime)/numbOfInterv;
    interpTime      = interpTime+timeStep;
    
    if (track_data(m+1,11)-track_data(m,11))>180,
        interpShearDir = track_data(m,11)+((track_data(m+1,11)-track_data(m,11))-360)*((1/timeStep)*interpTime)/numbOfInterv;
        if interpShearDir<0,
            interpShearDir = 360-interpShearDir;
        end
    elseif (track_data(m+1,11)-track_data(m,11))<-180,
        interpShearDir = track_data(m,11)+(360-(track_data(m+1,11)-track_data(m,11)))*((1/timeStep)*interpTime)/numbOfInterv;
        if interpShearDir>360,
            interpShearDir = interpShearDir-360;
        end
    else
        interpShearDir = track_data(m,11)+(track_data(m+1,11)-track_data(m,11))*((1/timeStep)*interpTime)/numbOfInterv;
    end
    
    
    interpShearDir =interpShearDir+180;
    if interpShearDir>360,
        interpShearDir=interpShearDir-360;
    end
    
    centerLon(index) = interpLon;
    centerLat(index) = interpLat;
    index=index+1;
    
% Definition of parameterization functions                      %
% TC Intensity based on best track wind speed (in kts)          %
% Shear amplitude from SHIPS (in kts)                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %interpShearAmpl=20.0;
    for r=1:50,
      rr=10*r-5;
      if interpWind<65.0,
        r_0(r) = 3.25*exp(-(rr-(interpRmax))^2/(125)^2)+0.8*exp(-(rr-(interpRmax+205))^2/(rr/2+25)^2)...
            -(3.25*exp(-(495-(interpRmax))^2/(125)^2)+0.8*exp(-(495-(interpRmax+205))^2/(495/2+25)^2));
        if interpShearAmpl<9.7,
          a_1(r) = -17.5*exp(-(rr-100)^2/((rr+50)*0.5)^2)+7.5;
          b_1(r) = 1.1*(rr-5)^0.6-exp((rr-125)/100);
          a_2(r) = 4-10*exp(-(rr-245)^2/250^2);
          b_2(r) = 1-10*exp(-(rr-335)^2/60^2);
        elseif ((interpShearAmpl>=9.7)&&(interpShearAmpl<19.4)),
          a_1(r) = -20*exp(-(rr-65)^2/((rr+50)*0.65)^2)+2.5;
          b_1(r) = 1.1*(rr-5)^0.75-((rr-25)/47.5)^2;
          a_2(r) = -2.5-12.5*exp(0-(rr-245)^2/100^2);
          b_2(r) = -1.5*(rr-5)^0.25;
        else
          a_1(r) = -35*exp(-(rr-85)^2/((rr+150)*0.35)^2)-5;
          b_1(r) = 1.1*(rr-5)^0.75-((rr-25)/50)^2+12.5;
          a_2(r) = -1.5-20*exp(-(rr-205)^2/100^2);
          b_2(r) = -2.5*(rr-5)^0.32;
        end
      elseif ((interpWind>=65.0)&&(interpWind<95.0)),
        r_0(r) = 7.0*exp(-(rr-interpRmax)^2/(100)^2)+1.5*exp(-(rr-(interpRmax+135))^2/(rr/2.25+25)^2)...
            -(7.0*exp(-(495-interpRmax)^2/(100)^2)+1.5*exp(-(495-(interpRmax+135))^2/(495/2.25+25)^2));
        if interpShearAmpl<9.7,
          a_1(r) = -15*exp(-(rr-100)^2/((rr+300)*0.25)^2)+(-8+rr/20)-exp((rr-175)/85);
          b_1(r) = 1.1*(rr-5)^0.35;
          a_2(r) = 10-8*exp((rr)^2/600^2);
          b_2(r) = 6.5-8*exp(-(rr-245)^2/250^2)-exp((rr-250)/100);
        elseif ((interpShearAmpl>=9.7)&&(interpShearAmpl<19.4)),
          a_1(r) = -15*exp(-(rr-100)^2/((rr+300)*0.2)^2)-9-exp((rr-250)/100);
          b_1(r) = 1.1*(rr-5)^0.75-2.5*(rr/50)^1.5+2.5;
          a_2(r) = -2.5-10*exp(-(rr-225)^2/(100-(rr-100)/10)^2);
          b_2(r) = -1.25*(rr-5)^0.2-exp((rr-400)/100);
        else
          a_1(r) = -25*exp(-(rr-100)^2/((rr+100)*0.75)^2)-9+exp((rr-200)/75);
          b_1(r) = 1.1*(rr-5)^0.7-((rr-25)/75)^2+5;
          a_2(r) = -1.5-16*exp(-(rr-375)^2/200^2);
          b_2(r) = -2.5*(rr-5)^0.23;
        end
      else
        r_0(r) = 12.5*exp(-(rr-interpRmax)^2/(75)^2)+3.25*exp(-(rr-(interpRmax+115))^2/(rr/3+25)^2)...
            -(12.5*exp(-(495-interpRmax)^2/(75)^2)+3.25*exp(-(495-(interpRmax+115))^2/(495/3+25)^2));
        if interpShearAmpl<9.7,
          a_1(r) = -15*exp(-(rr-100)^2/((rr+300)*0.5)^2)+(5+rr/20)-2*exp((rr-75)/150);
          b_1(r) = 1.1*(rr-5)^0.35;
          a_2(r) = 1.5-11.5*exp(-(rr-375)^2/125^2);
          b_2(r) = 6.5-8*exp(-(rr-245)^2/250^2)+22.5*exp(-(rr-410)^2/50^2);
        elseif ((interpShearAmpl>=9.7)&&(interpShearAmpl<19.4)),
          a_1(r) = -15*exp(-(rr-150)^2/((rr+500)*0.15)^2)-11;
          b_1(r) = 1.1*((rr-5)/10)^1.3-7*((rr-50)/100)^2;
          a_2(r) = 1-10*exp(-(rr-350)^2/125^2);
          b_2(r) = 4.5-20*exp(-(rr-300)^2/100^2);
        else
          a_1(r) = -15*exp(-(rr-100)^2/((rr+300)*0.35)^2)-15;
          b_1(r) = 1.1*(rr-5)^0.72-((rr-25)/50)^2+5;
          a_2(r) = -1.5-20*exp(-(rr-245)^2/50^2)+15*exp(-(rr-340)^2/50^2);
          b_2(r) = -5-10*exp(-(rr-245)^2/100^2);
        end
      end
    end

% Construction of rainfall footprint in cylindrical coordinates %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    for r=1:50,
      for theta=1:360,
        rain(r,theta)= r_0(r)+...
            a_1(r)*(r_0(r)/100)*cos(theta*2*pi/360)+b_1(r)*(r_0(r)/100)*sin(theta*2*pi/360)+...
            a_2(r)*(r_0(r)/100)*cos(2*theta*2*pi/360)+b_2(r)*(r_0(r)/100)*sin(2*theta*2*pi/360);
        %rain(r,theta)= a_1(r)*(r_0(r)/100)*cos(theta*2*pi/360)+b_1(r)*(r_0(r)/100)*sin(theta*2*pi/360)+...
        %   a_2(r)*(r_0(r)/100)*cos(2*theta*2*pi/360)+b_2(r)*(r_0(r)/100)*sin(2*theta*2*pi/360);
      end
    end
  
% footprint assigned lat/lon components                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    for r=1:50,
      rr=10*r-5;
      for theta=1:360,
        rainLon_north(r,theta) = rr*cos(theta*2*pi/360)/111.1214;
        rainLat_north(r,theta) = rr*sin(theta*2*pi/360)/111.1214;
      end
    end

    %interpShearDir=290;
    for r=1:50,
      for theta=1:360,
        rainLat(r,theta) = interpLat-rainLon_north(r,theta)*sin(interpShearDir*2*pi/360)+...
            rainLat_north(r,theta)*cos(interpShearDir*2*pi/360);
        rainLon(r,theta) = -interpLon+rainLon_north(r,theta)*cos(interpShearDir*2*pi/360)+...
            rainLat_north(r,theta)*sin(interpShearDir*2*pi/360);
      end
    end
    
% Interpolation of footprint to master geo-grid, a 1/10-degree resolution %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

interpRain=griddata(rainLon,rainLat,rain,masterLon,masterLat,'nearest');
master_history=master_history+interpRain/(1/timeStep);

% Accumulate the result into the master arrays, ultimately the output     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
  end
end

aek=master_history;
k=find(aek<0.1);%&aek>-5);
master_history(k)=NaN;


hf=figure;
set(hf,'paperorientation','portrait','paperposition',[1.0,1.0,9.0,6.5]);
%
ha1=axes('position',[.1,.1,.88,.88],'linewidth',2,'Fontsize',12,'fontname','Arial','box','off');
%axis equal;

hold on
whitebg('black');

v=[1,3,5,7,10,15,20,25];
%[c1,h1]=contourf(rainLon,rainLat,rain);
[c1,h1]=contourf(masterLon,masterLat,2*master_history/25.4,[50]);
set(h1,'edgecolor','none');
%set(h1,'linewidth',2);
%colormap hsv
%colormapeditor
hold on
coast=load('feasibility/coastline.txt');
hold on;
plot(coast(:,1),coast(:,2),'w');
hold on
for i=50:50:500,
  plot(masterLon(i,:), masterLat(i,:),'y--')
  hold on
end
for j=50:50:500,
  plot(masterLon(:,j), masterLat(:,j),'y--')
  hold on
end
hold on
plot(-centerLon,centerLat,'w.','markersize',3);
hold on

%axis([-95 -60, 20 50]);
%axis off
xlabel ('Longitude (degrees)');
ylabel ('Latitude (degrees)');

%caxis([1 5]);
gcf=colorbar(jet(256),'Fontsize',12,'fontname','Arial');
hold off;
