function [pdiland, jland, jlandhur, jlandmajor] = landfallpdi(lat, long, v, jmax)
%
% Determines power dissipation summed on first land fall points, using
% bathymetry data set. 
%    Copyright 2009 WindRiskTech L.L.C.
%    Last modified September, 2013
%-------------------------------------------------------------------------
load bathymetry.mat bathy
bathy=min(bathy,200); 
%
[n,~]=size(v);
jland=zeros(n,20);
jlandhur=zeros(n,20);
jlandmajor=zeros(n,20);
pdiland=zeros(1,n);
%
for i=1:n
    hold=-1000;
    jn=0;jnhur=0;jnmajor=0;
    for j=1:jmax(i)
        if long(i,j) >= 360
            long(i,j)=long(i,j)-360;
        end    
        ib=1+floor(4*long(i,j));
        if ib > 1440
            ib=1;
        end    
        ibp=ib+1;
        if ibp > 1440
            ibp=1;
        end    
        jb=1+floor(4*(lat(i,j)+90));
        b1=bathy(ib,jb);
        b2=bathy(ib,jb+1);
        b3=bathy(ibp,jb);
        b4=bathy(ibp,jb+1);
        dely=1+4*(lat(i,j)+90)-jb;
        delx=1+4*long(i,j)-ib;
        d1=(1.-delx)*(1.-dely);
        d2=dely*(1.-delx);
        d3=delx*(1.-dely);
        d4=delx*dely;
        h=-exp(d1.*log(-b1+220)+d2.*log(-b2+220)+d3.*log(-b3+220)+d4.*log(-b4+220))+220;
        %h=d1*b1+d2*b2+d3*b3+d4*b4;
        if hold < 0 && h > 0 
        %if hold < 0 && h > 0 && lat(i,j) >= 24
            jn=jn+1;
            jland(i,jn)=j;
            if v(i,j) >= 63
                jnhur=jnhur+1;
                jlandhur(i,jnhur)=j;
            end
            if v(i,j) >= 95
                jnmajor=jnmajor+1;
                jlandmajor(i,jnmajor)=j;
            end
            %
            jm=max(j-1,1);
            pdiland(i)=pdiland(i)+v(i,jm).^3;
            %pdiland(i)=pdiland(i)+damage(v(i,jm));
        end   
        hold=h;
    end
end    
