% Bree Aldridge, 2/15/2006
% compute DLEs using the trajectory end points in integration_results.mat

% load the integration results
load integration_results

% compute the 8-dimensional gradient 
% note that the gradient function in matlab switches the first two
% directions, so "A" and "B" corresponding to casp-8 and casp-8* are alway
% switched.  
[DAB,DAA,DAC,DAD,DAE,DAF,DAG,DAH]=gradient(stor(:,:,:,:,:,:,:,:,1),c8avar,c8var,c3var,c3avar,c8c3avar,c8ac3var,xiapvar,c3axvar);
[DBB,DBA,DBC,DBD,DBE,DBF,DBG,DBH]=gradient(stor(:,:,:,:,:,:,:,:,2),c8avar,c8var,c3var,c3avar,c8c3avar,c8ac3var,xiapvar,c3axvar);
[DCB,DCA,DCC,DCD,DCE,DCF,DCG,DCH]=gradient(stor(:,:,:,:,:,:,:,:,3),c8avar,c8var,c3var,c3avar,c8c3avar,c8ac3var,xiapvar,c3axvar);
[DDB,DDA,DDC,DDD,DDE,DDF,DDG,DDH]=gradient(stor(:,:,:,:,:,:,:,:,4),c8avar,c8var,c3var,c3avar,c8c3avar,c8ac3var,xiapvar,c3axvar);
[DEB,DEA,DEC,DED,DEE,DEF,DEG,DEH]=gradient(stor(:,:,:,:,:,:,:,:,5),c8avar,c8var,c3var,c3avar,c8c3avar,c8ac3var,xiapvar,c3axvar);
[DFB,DFA,DFC,DFD,DFE,DFF,DFG,DFH]=gradient(stor(:,:,:,:,:,:,:,:,6),c8avar,c8var,c3var,c3avar,c8c3avar,c8ac3var,xiapvar,c3axvar);
[DGB,DGA,DGC,DGD,DGE,DGF,DGG,DGH]=gradient(stor(:,:,:,:,:,:,:,:,7),c8avar,c8var,c3var,c3avar,c8c3avar,c8ac3var,xiapvar,c3axvar);
[DHB,DHA,DHC,DHD,DHE,DHF,DHG,DHH]=gradient(stor(:,:,:,:,:,:,:,:,8),c8avar,c8var,c3var,c3avar,c8c3avar,c8ac3var,xiapvar,c3axvar);

% open the file to write the space delimited results
fid=fopen('DLE_results','w+');
fprintf(fid,'%s\n',['c8 c8a c3 c3a c8c3a c8ac3 xiap c3axvar DLE']); %in the text file, annotate the initial consitions (phase-space location) and DLEs

% iterate through the grid points (now, the gradients) to compte DLEs and
% print the results
for a=1:length(c8var)
    for b=1:length(c8avar)       
        for c=1:length(c3var)           
            for d=1:length(c3avar)               
                for e=1:length(c8c3avar)                 
                    for f=1:length(c8ac3var)
                        for g=1:length(xiapvar)
                            for h=1:length(c3axvar) 
                                dxdxo=[DAA(a,b,c,d,e,f,g,h) DAB(a,b,c,d,e,f,g,h) DAC(a,b,c,d,e,f,g,h) DAD(a,b,c,d,e,f,g,h) DAE(a,b,c,d,e,f,g,h) DAF(a,b,c,d,e,f,g,h) DAG(a,b,c,d,e,f,g,h) DAH(a,b,c,d,e,f,g,h);
                                    DBA(a,b,c,d,e,f,g,h) DBB(a,b,c,d,e,f,g,h) DBC(a,b,c,d,e,f,g,h) DBD(a,b,c,d,e,f,g,h) DBE(a,b,c,d,e,f,g,h) DBF(a,b,c,d,e,f,g,h) DBG(a,b,c,d,e,f,g,h) DBH(a,b,c,d,e,f,g,h);
                                    DCA(a,b,c,d,e,f,g,h) DCB(a,b,c,d,e,f,g,h) DCC(a,b,c,d,e,f,g,h) DCD(a,b,c,d,e,f,g,h) DCE(a,b,c,d,e,f,g,h) DCF(a,b,c,d,e,f,g,h) DCG(a,b,c,d,e,f,g,h) DCH(a,b,c,d,e,f,g,h);
                                    DDA(a,b,c,d,e,f,g,h) DDB(a,b,c,d,e,f,g,h) DDC(a,b,c,d,e,f,g,h) DDD(a,b,c,d,e,f,g,h) DDE(a,b,c,d,e,f,g,h) DDF(a,b,c,d,e,f,g,h) DDG(a,b,c,d,e,f,g,h) DDH(a,b,c,d,e,f,g,h);
                                    DEA(a,b,c,d,e,f,g,h) DEB(a,b,c,d,e,f,g,h) DEC(a,b,c,d,e,f,g,h) DED(a,b,c,d,e,f,g,h) DEE(a,b,c,d,e,f,g,h) DEF(a,b,c,d,e,f,g,h) DEG(a,b,c,d,e,f,g,h) DEH(a,b,c,d,e,f,g,h);
                                    DFA(a,b,c,d,e,f,g,h) DFB(a,b,c,d,e,f,g,h) DFC(a,b,c,d,e,f,g,h) DFD(a,b,c,d,e,f,g,h) DFE(a,b,c,d,e,f,g,h) DFF(a,b,c,d,e,f,g,h) DFG(a,b,c,d,e,f,g,h) DFH(a,b,c,d,e,f,g,h);
                                    DGA(a,b,c,d,e,f,g,h) DGB(a,b,c,d,e,f,g,h) DGC(a,b,c,d,e,f,g,h) DGD(a,b,c,d,e,f,g,h) DGE(a,b,c,d,e,f,g,h) DGF(a,b,c,d,e,f,g,h) DGG(a,b,c,d,e,f,g,h) DGH(a,b,c,d,e,f,g,h);
                                    DHA(a,b,c,d,e,f,g,h) DHB(a,b,c,d,e,f,g,h) DHC(a,b,c,d,e,f,g,h) DHD(a,b,c,d,e,f,g,h) DHE(a,b,c,d,e,f,g,h) DHF(a,b,c,d,e,f,g,h) DHG(a,b,c,d,e,f,g,h) DHH(a,b,c,d,e,f,g,h)];
                                S=dxdxo'*dxdxo;
                                lmlam(a,b,c,d,e,f,g,h)=log(max(eig(S)));
                                fprintf(fid,'%s\n',[num2str(c8var(a)) ' ' num2str(c8avar(b)) ' ' num2str(c3var(c)) ' ' num2str(c3avar(d)) ' ' num2str(c8c3avar(e)) ' ' num2str(c8ac3var(f)) ' ' num2str(xiapvar(g)) ' ' num2str(c3axvar(h)) ' ' num2str(lmlam(a,b,c,d,e,f,g,h)) ' ']);
                            end
                        end
                    end
                end
            end
        end
    end
end 

fclose(fid);