text:
/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
/*                                                                  */
/* Linear Algebra Data and Routines File                            */
/*                                                                  */
/* Generated by KPP-2.2 symbolic chemistry Kinetics PreProcessor    */
/*       (http://www.cs.vt.edu/~asandu/Software/KPP)                */
/* KPP is distributed under GPL, the general public licence         */
/*       (http://www.gnu.org/copyleft/gpl.html)                     */
/* (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa           */
/* (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech            */
/*     With important contributions from:                           */
/*        M. Damian, Villanova University, USA                      */
/*        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany */
/*                                                                  */
/* File                 : saprc99_LinearAlgebra.c                   */
/* Time                 : Wed Jun 11 10:09:48 2008                  */
/* Working directory    : /Users/jlinford/workspace/saprc99         */
/* Equation file        : saprc99.kpp                               */
/* Output root filename : saprc99                                   */
/*                                                                  */
/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */

#define USE_SDK_BLAS    0

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "saprc99_Parameters.h"
#include "saprc99_Global.h"
#include "saprc99_Sparse.h"

#if DO_CHEMISTRY == 1

/* For SPU-optimized BLAS */
#if USE_SDK_BLAS == 1
#include <blas.h>
#endif

/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
/*                                                                  */
/* SPARSE_UTIL - SPARSE utility functions                           */
/*   Arguments :                                                    */
/*                                                                  */
/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */

int KppDecomp( double *JVS )
{
double W[74];
double a;
int k, kk, j, jj;

  for( k = 0; k < 74; k++ ) {
    if( JVS[ LU_DIAG[k] ] == 0.0 ) return k+1;
    for( kk = LU_CROW[k]; kk < LU_CROW[k+1]; kk++ )
      W[ LU_ICOL[kk] ] = JVS[kk];
    for( kk = LU_CROW[k]; kk < LU_DIAG[k]; kk++ ) {
      j = LU_ICOL[kk];
      a = -W[j] / JVS[ LU_DIAG[j] ];
      W[j] = -a;
      for( jj = LU_DIAG[j]+1; jj < LU_CROW[j+1]; jj++ )
        W[ LU_ICOL[jj] ] += a*JVS[jj];
    }
    for( kk = LU_CROW[k]; kk < LU_CROW[k+1]; kk++ )
      JVS[kk] = W[ LU_ICOL[kk] ];
  }
  return 0;
}

/* End of SPARSE_UTIL function                                      */
/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */

/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
/*                                                                  */
/* KppSolve - sparse back substitution                              */
/*   Arguments :                                                    */
/*      JVS       - sparse Jacobian of variables                    */
/*      X         - Vector for variables                            */
/*                                                                  */
/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */

void KppSolve(
  double JVS[],                           /* sparse Jacobian of variables */
  double X[]                              /* Vector for variables */
)
{
  X[21] = X[21]-JVS[91]*X[20];
  X[31] = X[31]-JVS[125]*X[23]-JVS[126]*X[30];
  X[32] = X[32]-JVS[129]*X[23]-JVS[130]*X[30];
  X[33] = X[33]-JVS[133]*X[23]-JVS[134]*X[30];
  X[34] = X[34]-JVS[138]*X[23]-JVS[139]*X[30];
  X[35] = X[35]-JVS[143]*X[27];
  X[37] = X[37]-JVS[154]*X[23]-JVS[155]*X[30];
  X[38] = X[38]-JVS[161]*X[30];
  X[39] = X[39]-JVS[167]*X[19]-JVS[168]*X[29]-JVS[169]*X[31]-JVS[170]
         *X[32]-JVS[171]*X[34];
  X[40] = X[40]-JVS[190]*X[23]-JVS[191]*X[30]-JVS[192]*X[31]-JVS[193]
         *X[32]-JVS[194]*X[33];
  X[41] = X[41]-JVS[202]*X[19]-JVS[203]*X[20]-JVS[204]*X[21]-JVS[205]
         *X[22]-JVS[206]*X[29];
  X[42] = X[42]-JVS[216]*X[17]-JVS[217]*X[33]-JVS[218]*X[35]-JVS[219]
         *X[37]-JVS[220]*X[38]-JVS[221]*X[40];
  X[44] = X[44]-JVS[242]*X[19]-JVS[243]*X[23]-JVS[244]*X[30]-JVS[245]
         *X[31]-JVS[246]*X[32]-JVS[247]*X[34]-JVS[248]*X[38]-JVS[249]
         *X[43];
  X[45] = X[45]-JVS[259]*X[33]-JVS[260]*X[38];
  X[47] = X[47]-JVS[276]*X[20]-JVS[277]*X[22]-JVS[278]*X[29]-JVS[279]
         *X[31]-JVS[280]*X[32]-JVS[281]*X[41]-JVS[282]*X[46];
  X[49] = X[49]-JVS[310]*X[46];
  X[51] = X[51]-JVS[322]*X[46];
  X[53] = X[53]-JVS[334]*X[46]-JVS[335]*X[52];
  X[54] = X[54]-JVS[341]*X[10]-JVS[342]*X[20]-JVS[343]*X[22]-JVS[344]
         *X[29]-JVS[345]*X[43]-JVS[346]*X[50]-JVS[347]*X[51]-JVS[348]
         *X[52];
  X[55] = X[55]-JVS[363]*X[19]-JVS[364]*X[20]-JVS[365]*X[22]-JVS[366]
         *X[25]-JVS[367]*X[26]-JVS[368]*X[28]-JVS[369]*X[29]-JVS[370]
         *X[41]-JVS[371]*X[43]-JVS[372]*X[44]-JVS[373]*X[45]-JVS[374]
         *X[46]-JVS[375]*X[48]-JVS[376]*X[49]-JVS[377]*X[50]-JVS[378]
         *X[51]-JVS[379]*X[52]-JVS[380]*X[53];
  X[56] = X[56]-JVS[399]*X[21]-JVS[400]*X[48]-JVS[401]*X[50]-JVS[402]
         *X[51]-JVS[403]*X[52];
  X[57] = X[57]-JVS[412]*X[9]-JVS[413]*X[43]-JVS[414]*X[46]-JVS[415]
         *X[48]-JVS[416]*X[49]-JVS[417]*X[50]-JVS[418]*X[52]-JVS[419]
         *X[53];
  X[58] = X[58]-JVS[426]*X[19]-JVS[427]*X[20]-JVS[428]*X[22]-JVS[429]
         *X[29]-JVS[430]*X[31]-JVS[431]*X[32]-JVS[432]*X[34]-JVS[433]
         *X[36]-JVS[434]*X[43]-JVS[435]*X[48]-JVS[436]*X[49]-JVS[437]
         *X[50]-JVS[438]*X[51]-JVS[439]*X[52]-JVS[440]*X[53]-JVS[441]
         *X[56]-JVS[442]*X[57];
  X[59] = X[59]-JVS[454]*X[20]-JVS[455]*X[22]-JVS[456]*X[29]-JVS[457]
         *X[49]-JVS[458]*X[50]-JVS[459]*X[51]-JVS[460]*X[52]-JVS[461]
         *X[53]-JVS[462]*X[56]-JVS[463]*X[57];
  X[60] = X[60]-JVS[474]*X[22]-JVS[475]*X[29]-JVS[476]*X[30]-JVS[477]
         *X[46]-JVS[478]*X[48]-JVS[479]*X[50]-JVS[480]*X[51]-JVS[481]
         *X[52]-JVS[482]*X[53]-JVS[483]*X[56]-JVS[484]*X[57];
  X[61] = X[61]-JVS[497]*X[34]-JVS[498]*X[43]-JVS[499]*X[46]-JVS[500]
         *X[48]-JVS[501]*X[49]-JVS[502]*X[50]-JVS[503]*X[51]-JVS[504]
         *X[52]-JVS[505]*X[53]-JVS[506]*X[57];
  X[62] = X[62]-JVS[517]*X[8]-JVS[518]*X[16]-JVS[519]*X[18]-JVS[520]
         *X[19]-JVS[521]*X[23]-JVS[522]*X[24]-JVS[523]*X[25]-JVS[524]
         *X[26]-JVS[525]*X[27]-JVS[526]*X[28]-JVS[527]*X[30]-JVS[528]
         *X[31]-JVS[529]*X[32]-JVS[530]*X[34]-JVS[531]*X[35]-JVS[532]
         *X[36]-JVS[533]*X[39]-JVS[534]*X[40]-JVS[535]*X[43]-JVS[536]
         *X[44]-JVS[537]*X[45]-JVS[538]*X[46]-JVS[539]*X[48]-JVS[540]
         *X[49]-JVS[541]*X[50]-JVS[542]*X[51]-JVS[543]*X[52]-JVS[544]
         *X[53]-JVS[545]*X[54]-JVS[546]*X[55]-JVS[547]*X[56]-JVS[548]
         *X[57]-JVS[549]*X[58]-JVS[550]*X[59]-JVS[551]*X[60]-JVS[552]
         *X[61];
  X[63] = X[63]-JVS[565]*X[19]-JVS[566]*X[20]-JVS[567]*X[22]-JVS[568]
         *X[23]-JVS[569]*X[29]-JVS[570]*X[30]-JVS[571]*X[46]-JVS[572]
         *X[48]-JVS[573]*X[50]-JVS[574]*X[51]-JVS[575]*X[52]-JVS[576]
         *X[53]-JVS[577]*X[56]-JVS[578]*X[57]-JVS[579]*X[58]-JVS[580]
         *X[59]-JVS[581]*X[60]-JVS[582]*X[61]-JVS[583]*X[62];
  X[64] = X[64]-JVS[595]*X[15]-JVS[596]*X[46]-JVS[597]*X[49]-JVS[598]
         *X[51]-JVS[599]*X[52]-JVS[600]*X[53]-JVS[601]*X[57]-JVS[602]
         *X[61]-JVS[603]*X[62]-JVS[604]*X[63];
  X[65] = X[65]-JVS[615]*X[21]-JVS[616]*X[25]-JVS[617]*X[29]-JVS[618]
         *X[41]-JVS[619]*X[43]-JVS[620]*X[46]-JVS[621]*X[48]-JVS[622]
         *X[50]-JVS[623]*X[52]-JVS[624]*X[53]-JVS[625]*X[54]-JVS[626]
         *X[56]-JVS[627]*X[57]-JVS[628]*X[58]-JVS[629]*X[59]-JVS[630]
         *X[60]-JVS[631]*X[61]-JVS[632]*X[62]-JVS[633]*X[63]-JVS[634]
         *X[64];
  X[66] = X[66]-JVS[644]*X[14]-JVS[645]*X[37]-JVS[646]*X[52]-JVS[647]
         *X[57]-JVS[648]*X[61]-JVS[649]*X[62]-JVS[650]*X[63]-JVS[651]
         *X[64]-JVS[652]*X[65];
  X[67] = X[67]-JVS[661]*X[10]-JVS[662]*X[19]-JVS[663]*X[20]-JVS[664]
         *X[22]-JVS[665]*X[23]-JVS[666]*X[29]-JVS[667]*X[30]-JVS[668]
         *X[31]-JVS[669]*X[32]-JVS[670]*X[33]-JVS[671]*X[34]-JVS[672]
         *X[36]-JVS[673]*X[38]-JVS[674]*X[43]-JVS[675]*X[45]-JVS[676]
         *X[46]-JVS[677]*X[48]-JVS[678]*X[49]-JVS[679]*X[50]-JVS[680]
         *X[51]-JVS[681]*X[52]-JVS[682]*X[53]-JVS[683]*X[56]-JVS[684]
         *X[57]-JVS[685]*X[58]-JVS[686]*X[59]-JVS[687]*X[60]-JVS[688]
         *X[61]-JVS[689]*X[62]-JVS[690]*X[63]-JVS[691]*X[64]-JVS[692]
         *X[65]-JVS[693]*X[66];
  X[68] = X[68]-JVS[701]*X[18]-JVS[702]*X[26]-JVS[703]*X[47]-JVS[704]
         *X[48]-JVS[705]*X[50]-JVS[706]*X[52]-JVS[707]*X[53]-JVS[708]
         *X[55]-JVS[709]*X[56]-JVS[710]*X[57]-JVS[711]*X[59]-JVS[712]
         *X[60]-JVS[713]*X[61]-JVS[714]*X[62]-JVS[715]*X[63]-JVS[716]
         *X[64]-JVS[717]*X[65]-JVS[718]*X[66]-JVS[719]*X[67];
  X[69] = X[69]-JVS[726]*X[12]-JVS[727]*X[13]-JVS[728]*X[14]-JVS[729]
         *X[15]-JVS[730]*X[17]-JVS[731]*X[18]-JVS[732]*X[21]-JVS[733]
         *X[24]-JVS[734]*X[26]-JVS[735]*X[27]-JVS[736]*X[35]-JVS[737]
         *X[42]-JVS[738]*X[44]-JVS[739]*X[45]-JVS[740]*X[46]-JVS[741]
         *X[47]-JVS[742]*X[48]-JVS[743]*X[49]-JVS[744]*X[50]-JVS[745]
         *X[51]-JVS[746]*X[52]-JVS[747]*X[53]-JVS[748]*X[54]-JVS[749]
         *X[55]-JVS[750]*X[56]-JVS[751]*X[57]-JVS[752]*X[58]-JVS[753]
         *X[59]-JVS[754]*X[60]-JVS[755]*X[61]-JVS[756]*X[62]-JVS[757]
         *X[63]-JVS[758]*X[64]-JVS[759]*X[65]-JVS[760]*X[66]-JVS[761]
         *X[67]-JVS[762]*X[68];
  X[70] = X[70]-JVS[768]*X[17]-JVS[769]*X[24]-JVS[770]*X[33]-JVS[771]
         *X[35]-JVS[772]*X[37]-JVS[773]*X[38]-JVS[774]*X[40]-JVS[775]
         *X[42]-JVS[776]*X[43]-JVS[777]*X[44]-JVS[778]*X[45]-JVS[779]
         *X[46]-JVS[780]*X[47]-JVS[781]*X[48]-JVS[782]*X[49]-JVS[783]
         *X[50]-JVS[784]*X[51]-JVS[785]*X[52]-JVS[786]*X[53]-JVS[787]
         *X[54]-JVS[788]*X[55]-JVS[789]*X[56]-JVS[790]*X[57]-JVS[791]
         *X[58]-JVS[792]*X[59]-JVS[793]*X[60]-JVS[794]*X[61]-JVS[795]
         *X[62]-JVS[796]*X[63]-JVS[797]*X[64]-JVS[798]*X[65]-JVS[799]
         *X[66]-JVS[800]*X[67]-JVS[801]*X[68]-JVS[802]*X[69];
  X[71] = X[71]-JVS[807]*X[11]-JVS[808]*X[12]-JVS[809]*X[23]-JVS[810]
         *X[29]-JVS[811]*X[31]-JVS[812]*X[32]-JVS[813]*X[40]-JVS[814]
         *X[41]-JVS[815]*X[48]-JVS[816]*X[49]-JVS[817]*X[50]-JVS[818]
         *X[51]-JVS[819]*X[52]-JVS[820]*X[53]-JVS[821]*X[54]-JVS[822]
         *X[56]-JVS[823]*X[57]-JVS[824]*X[58]-JVS[825]*X[59]-JVS[826]
         *X[60]-JVS[827]*X[61]-JVS[828]*X[62]-JVS[829]*X[63]-JVS[830]
         *X[64]-JVS[831]*X[65]-JVS[832]*X[66]-JVS[833]*X[67]-JVS[834]
         *X[68]-JVS[835]*X[69]-JVS[836]*X[70];
  X[72] = X[72]-JVS[840]*X[13]-JVS[841]*X[44]-JVS[842]*X[45]-JVS[843]
         *X[48]-JVS[844]*X[49]-JVS[845]*X[51]-JVS[846]*X[52]-JVS[847]
         *X[53]-JVS[848]*X[57]-JVS[849]*X[58]-JVS[850]*X[59]-JVS[851]
         *X[60]-JVS[852]*X[61]-JVS[853]*X[62]-JVS[854]*X[63]-JVS[855]
         *X[64]-JVS[856]*X[65]-JVS[857]*X[66]-JVS[858]*X[67]-JVS[859]
         *X[68]-JVS[860]*X[69]-JVS[861]*X[70]-JVS[862]*X[71];
  X[73] = X[73]-JVS[865]*X[8]-JVS[866]*X[9]-JVS[867]*X[10]-JVS[868]
         *X[16]-JVS[869]*X[18]-JVS[870]*X[19]-JVS[871]*X[20]-JVS[872]
         *X[22]-JVS[873]*X[23]-JVS[874]*X[24]-JVS[875]*X[25]-JVS[876]
         *X[28]-JVS[877]*X[29]-JVS[878]*X[30]-JVS[879]*X[31]-JVS[880]
         *X[32]-JVS[881]*X[33]-JVS[882]*X[34]-JVS[883]*X[36]-JVS[884]
         *X[37]-JVS[885]*X[38]-JVS[886]*X[39]-JVS[887]*X[40]-JVS[888]
         *X[41]-JVS[889]*X[42]-JVS[890]*X[43]-JVS[891]*X[44]-JVS[892]
         *X[45]-JVS[893]*X[46]-JVS[894]*X[48]-JVS[895]*X[49]-JVS[896]
         *X[50]-JVS[897]*X[51]-JVS[898]*X[52]-JVS[899]*X[53]-JVS[900]
         *X[54]-JVS[901]*X[55]-JVS[902]*X[56]-JVS[903]*X[57]-JVS[904]
         *X[58]-JVS[905]*X[59]-JVS[906]*X[60]-JVS[907]*X[61]-JVS[908]
         *X[62]-JVS[909]*X[63]-JVS[910]*X[64]-JVS[911]*X[65]-JVS[912]
         *X[66]-JVS[913]*X[67]-JVS[914]*X[68]-JVS[915]*X[69]-JVS[916]
         *X[70]-JVS[917]*X[71]-JVS[918]*X[72];
  X[73] = X[73]/JVS[919];
  X[72] = (X[72]-JVS[864]*X[73])/(JVS[863]);
  X[71] = (X[71]-JVS[838]*X[72]-JVS[839]*X[73])/(JVS[837]);
  X[70] = (X[70]-JVS[804]*X[71]-JVS[805]*X[72]-JVS[806]*X[73])
         /(JVS[803]);
  X[69] = (X[69]-JVS[764]*X[70]-JVS[765]*X[71]-JVS[766]*X[72]-JVS[767]
         *X[73])/(JVS[763]);
  X[68] = (X[68]-JVS[721]*X[69]-JVS[722]*X[70]-JVS[723]*X[71]-JVS[724]
         *X[72]-JVS[725]*X[73])/(JVS[720]);
  X[67] = (X[67]-JVS[695]*X[68]-JVS[696]*X[69]-JVS[697]*X[70]-JVS[698]
         *X[71]-JVS[699]*X[72]-JVS[700]*X[73])/(JVS[694]);
  X[66] = (X[66]-JVS[654]*X[67]-JVS[655]*X[68]-JVS[656]*X[69]-JVS[657]
         *X[70]-JVS[658]*X[71]-JVS[659]*X[72]-JVS[660]*X[73])
         /(JVS[653]);
  X[65] = (X[65]-JVS[636]*X[66]-JVS[637]*X[67]-JVS[638]*X[68]-JVS[639]
         *X[69]-JVS[640]*X[70]-JVS[641]*X[71]-JVS[642]*X[72]-JVS[643]
         *X[73])/(JVS[635]);
  X[64] = (X[64]-JVS[606]*X[65]-JVS[607]*X[66]-JVS[608]*X[67]-JVS[609]
         *X[68]-JVS[610]*X[69]-JVS[611]*X[70]-JVS[612]*X[71]-JVS[613]
         *X[72]-JVS[614]*X[73])/(JVS[605]);
  X[63] = (X[63]-JVS[585]*X[64]-JVS[586]*X[65]-JVS[587]*X[66]-JVS[588]
         *X[67]-JVS[589]*X[68]-JVS[590]*X[69]-JVS[591]*X[70]-JVS[592]
         *X[71]-JVS[593]*X[72]-JVS[594]*X[73])/(JVS[584]);
  X[62] = (X[62]-JVS[554]*X[63]-JVS[555]*X[64]-JVS[556]*X[65]-JVS[557]
         *X[66]-JVS[558]*X[67]-JVS[559]*X[68]-JVS[560]*X[69]-JVS[561]
         *X[70]-JVS[562]*X[71]-JVS[563]*X[72]-JVS[564]*X[73])
         /(JVS[553]);
  X[61] = (X[61]-JVS[508]*X[62]-JVS[509]*X[64]-JVS[510]*X[66]-JVS[511]
         *X[68]-JVS[512]*X[69]-JVS[513]*X[70]-JVS[514]*X[71]-JVS[515]
         *X[72]-JVS[516]*X[73])/(JVS[507]);
  X[60] = (X[60]-JVS[486]*X[61]-JVS[487]*X[63]-JVS[488]*X[65]-JVS[489]
         *X[66]-JVS[490]*X[67]-JVS[491]*X[68]-JVS[492]*X[69]-JVS[493]
         *X[70]-JVS[494]*X[71]-JVS[495]*X[72]-JVS[496]*X[73])
         /(JVS[485]);
  X[59] = (X[59]-JVS[465]*X[60]-JVS[466]*X[61]-JVS[467]*X[63]-JVS[468]
         *X[65]-JVS[469]*X[67]-JVS[470]*X[68]-JVS[471]*X[69]-JVS[472]
         *X[70]-JVS[473]*X[73])/(JVS[464]);
  X[58] = (X[58]-JVS[444]*X[59]-JVS[445]*X[60]-JVS[446]*X[61]-JVS[447]
         *X[62]-JVS[448]*X[63]-JVS[449]*X[67]-JVS[450]*X[68]-JVS[451]
         *X[69]-JVS[452]*X[70]-JVS[453]*X[73])/(JVS[443]);
  X[57] = (X[57]-JVS[421]*X[61]-JVS[422]*X[68]-JVS[423]*X[69]-JVS[424]
         *X[70]-JVS[425]*X[73])/(JVS[420]);
  X[56] = (X[56]-JVS[405]*X[57]-JVS[406]*X[61]-JVS[407]*X[63]-JVS[408]
         *X[68]-JVS[409]*X[69]-JVS[410]*X[70]-JVS[411]*X[73])
         /(JVS[404]);
  X[55] = (X[55]-JVS[382]*X[56]-JVS[383]*X[57]-JVS[384]*X[59]-JVS[385]
         *X[60]-JVS[386]*X[61]-JVS[387]*X[62]-JVS[388]*X[63]-JVS[389]
         *X[64]-JVS[390]*X[65]-JVS[391]*X[66]-JVS[392]*X[67]-JVS[393]
         *X[68]-JVS[394]*X[69]-JVS[395]*X[70]-JVS[396]*X[71]-JVS[397]
         *X[72]-JVS[398]*X[73])/(JVS[381]);
  X[54] = (X[54]-JVS[350]*X[56]-JVS[351]*X[57]-JVS[352]*X[58]-JVS[353]
         *X[59]-JVS[354]*X[60]-JVS[355]*X[61]-JVS[356]*X[64]-JVS[357]
         *X[66]-JVS[358]*X[68]-JVS[359]*X[70]-JVS[360]*X[71]-JVS[361]
         *X[72]-JVS[362]*X[73])/(JVS[349]);
  X[53] = (X[53]-JVS[337]*X[57]-JVS[338]*X[61]-JVS[339]*X[70]-JVS[340]
         *X[73])/(JVS[336]);
  X[52] = (X[52]-JVS[330]*X[57]-JVS[331]*X[61]-JVS[332]*X[70]-JVS[333]
         *X[73])/(JVS[329]);
  X[51] = (X[51]-JVS[324]*X[52]-JVS[325]*X[57]-JVS[326]*X[61]-JVS[327]
         *X[70]-JVS[328]*X[73])/(JVS[323]);
  X[50] = (X[50]-JVS[318]*X[57]-JVS[319]*X[61]-JVS[320]*X[70]-JVS[321]
         *X[73])/(JVS[317]);
  X[49] = (X[49]-JVS[312]*X[52]-JVS[313]*X[57]-JVS[314]*X[61]-JVS[315]
         *X[70]-JVS[316]*X[73])/(JVS[311]);
  X[48] = (X[48]-JVS[306]*X[57]-JVS[307]*X[61]-JVS[308]*X[70]-JVS[309]
         *X[73])/(JVS[305]);
  X[47] = (X[47]-JVS[284]*X[48]-JVS[285]*X[50]-JVS[286]*X[52]-JVS[287]
         *X[53]-JVS[288]*X[56]-JVS[289]*X[57]-JVS[290]*X[59]-JVS[291]
         *X[60]-JVS[292]*X[61]-JVS[293]*X[62]-JVS[294]*X[63]-JVS[295]
         *X[64]-JVS[296]*X[65]-JVS[297]*X[66]-JVS[298]*X[67]-JVS[299]
         *X[68]-JVS[300]*X[69]-JVS[301]*X[70]-JVS[302]*X[71]-JVS[303]
         *X[72]-JVS[304]*X[73])/(JVS[283]);
  X[46] = (X[46]-JVS[272]*X[57]-JVS[273]*X[61]-JVS[274]*X[70]-JVS[275]
         *X[73])/(JVS[271]);
  X[45] = (X[45]-JVS[262]*X[62]-JVS[263]*X[64]-JVS[264]*X[66]-JVS[265]
         *X[68]-JVS[266]*X[69]-JVS[267]*X[70]-JVS[268]*X[71]-JVS[269]
         *X[72]-JVS[270]*X[73])/(JVS[261]);
  X[44] = (X[44]-JVS[251]*X[45]-JVS[252]*X[48]-JVS[253]*X[51]-JVS[254]
         *X[57]-JVS[255]*X[61]-JVS[256]*X[62]-JVS[257]*X[70]-JVS[258]
         *X[73])/(JVS[250]);
  X[43] = (X[43]-JVS[238]*X[57]-JVS[239]*X[61]-JVS[240]*X[70]-JVS[241]
         *X[73])/(JVS[237]);
  X[42] = (X[42]-JVS[223]*X[44]-JVS[224]*X[45]-JVS[225]*X[49]-JVS[226]
         *X[51]-JVS[227]*X[52]-JVS[228]*X[53]-JVS[229]*X[54]-JVS[230]
         *X[55]-JVS[231]*X[58]-JVS[232]*X[61]-JVS[233]*X[62]-JVS[234]
         *X[69]-JVS[235]*X[70]-JVS[236]*X[73])/(JVS[222]);
  X[41] = (X[41]-JVS[208]*X[48]-JVS[209]*X[50]-JVS[210]*X[52]-JVS[211]
         *X[56]-JVS[212]*X[61]-JVS[213]*X[69]-JVS[214]*X[70]-JVS[215]
         *X[73])/(JVS[207]);
  X[40] = (X[40]-JVS[196]*X[49]-JVS[197]*X[51]-JVS[198]*X[53]-JVS[199]
         *X[61]-JVS[200]*X[70]-JVS[201]*X[73])/(JVS[195]);
  X[39] = (X[39]-JVS[173]*X[40]-JVS[174]*X[43]-JVS[175]*X[44]-JVS[176]
         *X[46]-JVS[177]*X[48]-JVS[178]*X[49]-JVS[179]*X[50]-JVS[180]
         *X[51]-JVS[181]*X[52]-JVS[182]*X[53]-JVS[183]*X[54]-JVS[184]
         *X[55]-JVS[185]*X[57]-JVS[186]*X[58]-JVS[187]*X[61]-JVS[188]
         *X[70]-JVS[189]*X[73])/(JVS[172]);
  X[38] = (X[38]-JVS[163]*X[45]-JVS[164]*X[62]-JVS[165]*X[70]-JVS[166]
         *X[73])/(JVS[162]);
  X[37] = (X[37]-JVS[157]*X[52]-JVS[158]*X[61]-JVS[159]*X[70]-JVS[160]
         *X[73])/(JVS[156]);
  X[36] = (X[36]-JVS[150]*X[62]-JVS[151]*X[63]-JVS[152]*X[67]-JVS[153]
         *X[73])/(JVS[149]);
  X[35] = (X[35]-JVS[145]*X[45]-JVS[146]*X[62]-JVS[147]*X[69]-JVS[148]
         *X[70])/(JVS[144]);
  X[34] = (X[34]-JVS[141]*X[61]-JVS[142]*X[73])/(JVS[140]);
  X[33] = (X[33]-JVS[136]*X[70]-JVS[137]*X[73])/(JVS[135]);
  X[32] = (X[32]-JVS[132]*X[73])/(JVS[131]);
  X[31] = (X[31]-JVS[128]*X[73])/(JVS[127]);
  X[30] = (X[30]-JVS[124]*X[73])/(JVS[123]);
  X[29] = (X[29]-JVS[122]*X[73])/(JVS[121]);
  X[28] = (X[28]-JVS[117]*X[63]-JVS[118]*X[65]-JVS[119]*X[67]-JVS[120]
         *X[73])/(JVS[116]);
  X[27] = (X[27]-JVS[112]*X[35]-JVS[113]*X[62]-JVS[114]*X[69]-JVS[115]
         *X[70])/(JVS[111]);
  X[26] = (X[26]-JVS[108]*X[55]-JVS[109]*X[62]-JVS[110]*X[68])
         /(JVS[107]);
  X[25] = (X[25]-JVS[104]*X[62]-JVS[105]*X[65]-JVS[106]*X[73])
         /(JVS[103]);
  X[24] = (X[24]-JVS[100]*X[62]-JVS[101]*X[69]-JVS[102]*X[73])
         /(JVS[99]);
  X[23] = (X[23]-JVS[98]*X[73])/(JVS[97]);
  X[22] = (X[22]-JVS[96]*X[73])/(JVS[95]);
  X[21] = (X[21]-JVS[93]*X[69]-JVS[94]*X[73])/(JVS[92]);
  X[20] = (X[20]-JVS[90]*X[73])/(JVS[89]);
  X[19] = (X[19]-JVS[88]*X[73])/(JVS[87]);
  X[18] = (X[18]-JVS[85]*X[68]-JVS[86]*X[73])/(JVS[84]);
  X[17] = (X[17]-JVS[82]*X[69]-JVS[83]*X[70])/(JVS[81]);
  X[16] = (X[16]-JVS[79]*X[62]-JVS[80]*X[73])/(JVS[78]);
  X[15] = (X[15]-JVS[76]*X[64]-JVS[77]*X[69])/(JVS[75]);
  X[14] = (X[14]-JVS[73]*X[66]-JVS[74]*X[69])/(JVS[72]);
  X[13] = (X[13]-JVS[70]*X[69]-JVS[71]*X[72])/(JVS[69]);
  X[12] = (X[12]-JVS[67]*X[69]-JVS[68]*X[71])/(JVS[66]);
  X[11] = (X[11]-JVS[62]*X[23]-JVS[63]*X[48]-JVS[64]*X[61]-JVS[65]
         *X[73])/(JVS[61]);
  X[10] = (X[10]-JVS[60]*X[73])/(JVS[59]);
  X[9] = (X[9]-JVS[58]*X[61])/(JVS[57]);
  X[8] = (X[8]-JVS[56]*X[73])/(JVS[55]);
  X[7] = (X[7]-JVS[52]*X[27]-JVS[53]*X[37]-JVS[54]*X[69])/(JVS[51]);
  X[6] = (X[6]-JVS[49]*X[27]-JVS[50]*X[69])/(JVS[48]);
  X[5] = (X[5]-JVS[44]*X[62]-JVS[45]*X[64]-JVS[46]*X[66]-JVS[47]*X[72])
        /(JVS[43]);
  X[4] = (X[4]-JVS[41]*X[62]-JVS[42]*X[71])/(JVS[40]);
  X[3] = (X[3]-JVS[27]*X[46]-JVS[28]*X[48]-JVS[29]*X[50]-JVS[30]*X[51]
        -JVS[31]*X[52]-JVS[32]*X[61]-JVS[33]*X[62]-JVS[34]*X[63]
        -JVS[35]*X[64]-JVS[36]*X[65]-JVS[37]*X[66]-JVS[38]*X[67]
        -JVS[39]*X[72])/(JVS[26]);
  X[2] = (X[2]-JVS[18]*X[50]-JVS[19]*X[52]-JVS[20]*X[61]-JVS[21]*X[62]
        -JVS[22]*X[63]-JVS[23]*X[65]-JVS[24]*X[67]-JVS[25]*X[71])
        /(JVS[17]);
  X[1] = (X[1]-JVS[4]*X[19]-JVS[5]*X[26]-JVS[6]*X[43]-JVS[7]*X[46]
        -JVS[8]*X[48]-JVS[9]*X[49]-JVS[10]*X[50]-JVS[11]*X[51]-JVS[12]
        *X[52]-JVS[13]*X[53]-JVS[14]*X[61]-JVS[15]*X[68]-JVS[16]*X[73])
        /(JVS[3]);
  X[0] = (X[0]-JVS[1]*X[8]-JVS[2]*X[73])/(JVS[0]);
}

/* End of KppSolve function                                         */
/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */

/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
/*                                                                  */
/* BLAS_UTIL - BLAS-LIKE utility functions                          */
/*   Arguments :                                                    */
/*                                                                  */
/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */

/*--------------------------------------------------------------

  BLAS/LAPACK-like subroutines used by the integration algorithms
  It is recommended to replace them by calls to the optimized
       BLAS/LAPACK library for your machine

   (C) Adrian Sandu, Aug. 2004

--------------------------------------------------------------*/

#define ZERO  (double)0.0
#define ONE   (double)1.0
#define HALF  (double)0.5
#define TWO   (double)2.0
#define MOD(A,B) (int)((A)%(B))

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
void WCOPY(int N, double X[], int incX, volatile double Y[], int incY)
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    copies a vector, x, to a vector, y:  y <- x
    only for incX=incY=1
    after BLAS
    replace this by the function from the optimized BLAS implementation:
        CALL  SCOPY(N,X,1,Y,1)   or   CALL  DCOPY(N,X,1,Y,1)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
{
#if USE_SDK_BLAS == 1
    dcopy_spu(X, Y, N);
#else
      int i, M;
      if (N <= 0) return;

      M = MOD(N,8);
      if( M != 0 ) {
        for ( i = 0; i < M; i++ )
          Y[i] = X[i];
        if( N < 8 ) return;
      }
      for ( i = M; i<N; i+=8 ) {
        Y[i] = X[i];
        Y[i + 1] = X[i + 1];
        Y[i + 2] = X[i + 2];
        Y[i + 3] = X[i + 3];
        Y[i + 4] = X[i + 4];
        Y[i + 5] = X[i + 5];
        Y[i + 6] = X[i + 6];
        Y[i + 7] = X[i + 7];
      }
#endif
} /* end function WCOPY */

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
void WAXPY(int N, double Alpha, double X[], int incX, double Y[], int incY )
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    constant times a vector plus a vector: y <- y + Alpha*x
    only for incX=incY=1
    after BLAS
    replace this by the function from the optimized BLAS implementation:
        CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
{
#if USE_SDK_BLAS == 1
    daxpy_spu(X, Y, Alpha, N);
#else
      int i, M;

      if (Alpha == ZERO) return;
      if (N  <=  0) return;

      M = MOD(N,4);
      if( M != 0 ) {
        for ( i = 0; i < M; i++ )
          Y[i] = Y[i] + Alpha*X[i];
        if ( N < 4 ) return;
      }

      for ( i = M; i < N; i += 4 ) {
        Y[i] = Y[i] + Alpha*X[i];
        Y[i + 1] = Y[i + 1] + Alpha*X[i + 1];
        Y[i + 2] = Y[i + 2] + Alpha*X[i + 2];
        Y[i + 3] = Y[i + 3] + Alpha*X[i + 3];
      }
#endif

} /* end function  WAXPY */

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
void WSCAL(int N, double Alpha, double X[], int incX)
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    constant times a vector: x(1:N) <- Alpha*x(1:N)
    only for incX=incY=1
    after BLAS
    replace this by the function from the optimized BLAS implementation:
        CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
{
#if USE_SDK_BLAS == 1
    dscal_spu(X, Alpha, N);
#else
      int i, M;

      if (Alpha == ONE) return;
      if (N  <=  0) return;

      M = MOD(N,5);
      if( M  !=  0 ) {
        if (Alpha == (-ONE))
          for ( i = 0; i < M; i++ )  X[i] = -X[i];
        else {
	  if (Alpha == ZERO)
            for ( i = 0; i < M; i++ ) X[i] = ZERO;
          else
            for ( i = 0; i < M; i++ ) X[i] = Alpha*X[i];
        } // end else
        if( N < 5 ) return;
      } // end if

      if (Alpha == (-ONE))
        for ( i = M; i<N; i+=5 ) {
          X[i]     = -X[i];
          X[i + 1] = -X[i + 1];
          X[i + 2] = -X[i + 2];
          X[i + 3] = -X[i + 3];
          X[i + 4] = -X[i + 4];
        } // end for
      else {
        if (Alpha == ZERO)
          for ( i = M; i < N; i += 5 ) {
            X[i]     = ZERO;
            X[i + 1] = ZERO;
            X[i + 2] = ZERO;
            X[i + 3] = ZERO;
            X[i + 4] = ZERO;
          } // end for
        else
          for ( i = M; i < N; i += 5 ) {
            X[i]     = Alpha*X[i];
            X[i + 1] = Alpha*X[i + 1];
            X[i + 2] = Alpha*X[i + 2];
            X[i + 3] = Alpha*X[i + 3];
            X[i + 4] = Alpha*X[i + 4];
           } // end for
      }  // else
#endif

} /* end function WSCAL */

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
double WLAMCH( char C )
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    returns epsilon machine
    after LAPACK
    replace this by the function from the optimized LAPACK implementation:
         CALL SLAMCH('E') or CALL DLAMCH('E')
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
{
      int i;
      double Suma;
      static double Eps;
      static char First = 1;

      if (First) {
        First = 0;
        Eps = pow(HALF,16);
        for ( i = 17; i <= 80; i++ ) {
          Eps = Eps*HALF;
	  Suma = ONE + Eps;
	  if (Suma <= ONE) break;
        } /* end for */
        if (i==80) {
	   printf("\nERROR IN WLAMCH. Very small EPS = %g\n",Eps);
           return (double)2.2e-16;
	}
        Eps *= TWO; i--;
      } /* end if First */

      return Eps;

} /* end function WLAMCH */

#endif

max_stars_repo_name: jlinford/fixedgrid

max_issues_repo_name: jlinford/fixedgrid

max_forks_repo_name: jlinford/fixedgrid


text:
/*******************************************************************************
NAME: diffimage

ALGORITHM DESCRIPTION:
  Calculates statistics in each input file
  Calculates PSNR between the two input files
  Performs an fftMatch between the two input files
  Checks for differences in stats, a PSNR that is too low, and for
    shifts in geolocation

ISSUES:
        The images must have at least 75% overlap.
        diffimage will occasionally fail to find the
        correct offset especially with noisy images
        or images with flat topography (low non-noise
        spatial frequency.)

*******************************************************************************/
#include "asf.h"
#include "asf_nan.h"
#include <unistd.h>
#include <math.h>
#include <ctype.h>
#include "fft.h"
#include "fft2d.h"
#include "ddr.h"
#include "asf_raster.h"
#include "typlim.h"
#include "float_image.h"
#include "uint8_image.h"
#include "proj.h"
#include "libasf_proj.h"
#include "asf_jpeg.h"
#include "asf_tiff.h"
#include "geo_tiffp.h"
#include "geo_keyp.h"
#include <png.h>
#include <gsl/gsl_math.h>
#include "diffimage_tolerances.h"
#include "geotiff_support.h"
#include "asf_complex.h"

#ifndef png_jmpbuf
#  define png_jmpbuf (png_ptr)    ((png_ptr)->jmpbuf)
#endif

#define MISSING_PSNR -32000
#define FILE1_FFTFILE "tmp_file1.img"
#define FILE1_FFTFILE_META "tmp_file1.meta"
#define FILE2_FFTFILE "tmp_file2.img"
#define FILE2_FFTFILE_META "tmp_file2.meta"
#define CORR_FILE "tmp_corr"

/**** PROTOTYPES ****/
void usage(char *name);

int main(int argc, char **argv)
{
  char *inFile1,*inFile2;
  char *outputFile = NULL;
  extern int optind;            /* argv index of the next argument */
  extern char *optarg;          /* current argv[] */
  int c;                        /* option letter from getopt() */
  char ch;
  extern FILE *fLog;            /* output file descriptor, stdout or log file */
  FILE *fError;                 /* Error log, stderr or log file */
  extern int logflag, quietflag;
  int outputflag=0;
  int bandflag=0;
  int strictflag=0;
  int band=0;
  char msg[1024];

  fLog = NULL;
  logflag = quietflag=0;

  /* process command line */
  // FIXME: Might want to add -band1 and -band2 flags that would allow comparing
  // any arbitrary band in file1 to any arbitrary band in file2
  while ((c=getopt(argc,argv,"o:l:b:s:")) != EOF)
  {
    ch = (char)c;
    switch (ch) {
      case 'l':/* -log <filename>, get logfile; this is sorta hacked */
        if (0==strncmp(optarg,"og",2)) {
          sscanf(argv[optind++], "%s", logFile);
          logflag=1;
          fLog = FOPEN(logFile, "w");
        }
        else {
          FREE(outputFile);
          usage(argv[0]);
        }
        break;
      case 'o':/* -output <filename>, get output filename ...
		  empty if no differences */
        if (0==strncmp(optarg,"utput",5)) {
	  outputFile= (char*) CALLOC(1024, sizeof(char));
	  sscanf(argv[optind++], "%s", outputFile);
          outputflag=1;
        }
        else {
          strcpy(outputFile, "");
          outputflag=0;
        }
        break;
      case 'b':/* -band flag */
        if (0==strncmp(optarg,"and",3)) {
          sscanf(argv[optind++], "%d", &band);
          bandflag=1;
        }
        else {
          FREE(outputFile);
          usage(argv[0]);
        }
        break;
      case 's': /* -strict flag */
        if (0==strncmp(optarg,"trict",5)) {
          strictflag=1;
        }
        else {
          FREE(outputFile);
          usage(argv[0]);
        }
        break;
      default:
        FREE(outputFile);
        usage(argv[0]);
        break;
    }
  }
  asfSplashScreen(argc, argv);

  // After parsing out the command line arguments, there should be 2 arguments
  // left... the two files to compare
  if ((argc-optind) != 2) {
    if ((argc-optind) > 2) {
      printf("\n => Too many inputs.\n");
    }
    if ((argc-optind) < 2) {
      printf("\n => Too few inputs.\n");
    }
    FREE(outputFile);
    usage(argv[0]);
  }
  else {
    // Grab the file names
    inFile1=argv[optind];
    inFile2=argv[optind+1];
  }
  if (strcmp(inFile1, inFile2) == 0) {
    asfPrintError("inFile1 and inFile2 must be different files\n");
  }
  // Set up output redirection for error and log messages
  if (logflag == 0) {
    fLog = NULL;
    fError = NULL;
  }
  else {
    fError = fLog;
  }
  if (!outputflag) {
      sprintf(msg, "Missing output file name ...file differences will be "
	      "directed to stderr (only)\n");
      printf("** Warning: ********\n%s** End of warning **\n\n", msg);
  }
  if (outputflag && strcmp(logFile, outputFile) == 0) {
    sprintf(msg, "Log file cannot be the same as the output file:\n"
            "     Log file: %s\n  Output file: %s\n",
            logFile, outputFile);
    if (outputFile) FREE(outputFile);
    asfPrintError(msg);
  }

  char **bands1 = NULL, **bands2 = NULL;
  int num_bands1, num_bands2, complex = 0;
  stats_t *stats1 = NULL, *stats2 = NULL;
  complex_stats_t *complex_stats1 = NULL, *complex_stats2 = NULL;
  psnr_t *psnrs = NULL;
  complex_psnr_t *complex_psnr = NULL;
  shift_data_t *data_shift = NULL;
  diffimage(inFile1, inFile2, outputFile, logFile, &bands1, &bands2,
	    &num_bands1, &num_bands2, &complex, &stats1, &stats2,
	    &complex_stats1, &complex_stats2, &psnrs, &complex_psnr,
	    &data_shift);

  return (0);
}

void usage(char *name)
{
  printf("\nUSAGE:\n"
      "   %s [-output <diff_output_file>] [-log <file>] [-band <band number>]\n"
      "             [-strict] <img1.ext> <img2.ext>\n"
         "\nOPTIONS:\n"
      "   -output <diff_output_file>:  output to write image differencing\n"
      "                 results to (required.)\n"
      "   -log <file>:  allows the output to be written to a log file\n"
      "                 in addition to stdout (not required but strongly suggested.)\n"
      "   -band <band number>:  allows comparing one selected band from within\n"
      "                 multi-banded image files, including the separate color bands\n"
      "                 found within color images.  Default is to compare all\n"
      "                 available bands.\n"
      "   -strict:      forces statistics comparison to utilize all available\n"
      "                 fields (min, max, mean, sdev, PSNR).  Without the -strict\n"
      "                 option specified, only mean, sdev, and PSNR are utilized.\n"
      "                 default is non-strict.\n"
      "\nINPUTS:\n"
      "   <img1.ext>:   any supported single-banded graphics file image.\n"
      "                 supported types include TIFF, GEOTIFF, PNG, PGM, PPM,\n"
      "                 JPEG, and ASF IMG.  File extension is REQUIRED.\n"
      "   <img2.ext>:   an image to line up with the first.\n"
      "                 image 2 should be the same graphics file type and not\n"
      "                 be bigger than image 1.  File extension is REQUIRED.\n"
      "\nDESCRIPTION:\n"
      "   1. diffimage calculates image statistics within each input image\n"
      "      and calculates the peak signal-to-noise (PSNR) between the two\n"
      "      images.\n"
      "   2. diffimage then lines up the two images, to slightly better\n"
      "      than single-pixel precision, and determines if any geolocation\n"
      "      shift has occurred between the two and the size of the shift.\n"
      "      Because an fft-based match is utilized it will work with images of\n"
      "      any size, but is most efficient when the image dimensions are\n"
      "      near a power of 2.  The images need not be square.\n"
      "   3. diffimage then compares image statistics and geolocation shift\n"
      "      (if it occurred) and determines if the two images are different from\n"
      "      each other or not.\n"
      "   4. If there are no differences, the output file will exist but will be\n"
      "      of zero length.  Otherwise, a summary of the differences will be placed\n"
      "      in both the output file and the log file (if specified.)\n\n"
          "Version:\n   " SVN_REV " (part of " TOOL_SUITE_NAME " "MAPREADY_VERSION_STRING ")\n\n",
      name);
  exit(1);
}

ext: c

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
/* dwt_undec.h
 * Paul Demorest, 2007/10
 */
#ifndef _DWT_UNDEC_H_
#define _DWT_UNDEC_H_
#include <gsl/gsl_wavelet.h>
#ifdef __cplusplus
extern "C" {
#endif
int dwt_undec_transform(double *in, double *out, size_t n, const gsl_wavelet *w);
int dwt_undec_inverse(double *in, double *out, size_t n, const gsl_wavelet *w);
#ifdef __cplusplus
}
#endif
#endif

ext: h


text:
/**
 *
 * @file core_csetvar.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mark Gates
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:49 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex32_t
 *
 * CORE_csetvar sets a single variable, x := alpha.
 *
 *******************************************************************************
 *
 * @param[in] alpha
 *         Scalar to set x to, passed by pointer so it can depend on runtime value.
 *
 * @param[out] x
 *         On exit, x = alpha.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_csetvar = PCORE_csetvar
#define CORE_csetvar PCORE_csetvar
#endif
void CORE_csetvar(const PLASMA_Complex32_t *alpha, PLASMA_Complex32_t *x)
{
    *x = *alpha;
}


text:
/*****************************************************************\
           __
          / /
		 / /                     __  __
		/ /______    _______    / / / / ________   __       __
	   / ______  \  /_____  \  / / / / / _____  | / /      / /
	  / /      | / _______| / / / / / / /____/ / / /      / /
	 / /      / / / _____  / / / / / / _______/ / /      / /
	/ /      / / / /____/ / / / / / / |______  / |______/ /
   /_/      /_/ |________/ / / / /  \_______/  \_______  /
                          /_/ /_/                     / /
			                                         / /
		       High Level Game Framework            /_/

  ---------------------------------------------------------------

  Copyright (c) 2007-2011 - Rodrigo Braz Monteiro.
  This file is subject to the terms of halley_license.txt.

\*****************************************************************/

#pragma once

#include <cmath>
#include <iostream>
#include "angle.h"
#include <halley/utils/utils.h>
#include "halley/text/string_converter.h"
#include <gsl/gsl>

#ifdef min
#undef min
#endif
#ifdef max
#undef max
#endif

namespace Halley {
	//////////////////////////////
	// Vector2D class declaration
	template <typename T=float, class U=Angle<float>>
	class Vector2D {
	private:
		constexpr T mod(T a, T b) const { return a % b; }

	public:
		T x,y;

		// Constructors
		constexpr Vector2D () : x(0), y(0) {}

		constexpr Vector2D (T _x, T _y) : x(_x), y(_y) {}

		template <typename V>
		constexpr explicit Vector2D (Vector2D<V> vec) : x(T(vec.x)), y(T(vec.y)) {}

		Vector2D (T length, U angle)
		{
			float s, c;
			angle.sincos(s, c);
			x = s*length;
			y = c*length;
		}

		// Getter
		constexpr inline T& operator[](size_t n)
		{
			Expects(n <= 1);
			return (&x)[n];
		}
		constexpr inline T operator[](size_t n) const
		{
			Expects(n <= 1);
			return (&x)[n];
		}

		// Assignment and comparison
		constexpr inline Vector2D& operator = (Vector2D param) { x = param.x; y = param.y; return *this; }
		constexpr inline Vector2D& operator = (T param) { x = param; y = param; return *this; }
		constexpr inline bool operator == (Vector2D param) const { return x == param.x && y == param.y; }
		constexpr inline bool operator != (Vector2D param) const { return x != param.x || y != param.y; }
		constexpr inline bool operator < (Vector2D param) const { return y != param.y ? y < param.y : x < param.x; }

		// Basic algebra
		constexpr inline Vector2D operator + (Vector2D param) const { return Vector2D(x + param.x,y + param.y); }
		constexpr inline Vector2D operator - (Vector2D param) const { return Vector2D(x - param.x,y - param.y); }
		constexpr inline Vector2D operator * (Vector2D param) const { return Vector2D(x * param.x,y * param.y); }
		constexpr inline Vector2D operator / (Vector2D param) const { return Vector2D(x / param.x,y / param.y); }
		constexpr inline Vector2D operator % (Vector2D param) const { return Vector2D(mod(x, param.x), mod(y, param.y)); }

		constexpr inline Vector2D modulo(Vector2D param) const { return Vector2D(Halley::modulo<T>(x, param.x), Halley::modulo<T>(y, param.y)); }
		constexpr inline Vector2D floorDiv(Vector2D param) const { return Vector2D(Halley::floorDiv(x, param.x), Halley::floorDiv(y, param.y)); }

		constexpr inline Vector2D operator - () const { return Vector2D(-x,-y); }

		template <typename V>
		constexpr inline Vector2D operator * (V param) const { return Vector2D(T(x * param), T(y * param)); }

		template <typename V>
		constexpr inline Vector2D operator / (V param) const { return Vector2D(T(x / param), T(y / param)); }

		// In-place operations
		constexpr inline Vector2D& operator += (Vector2D param) { x += param.x; y += param.y; return *this; }
		constexpr inline Vector2D& operator -= (Vector2D param) { x -= param.x; y -= param.y; return *this; }
		constexpr inline Vector2D& operator *= (const T param) { x *= param; y *= param; return *this; }
		constexpr inline Vector2D& operator /= (const T param) { x /= param; y /= param; return *this; }

		// Get the normalized vector (unit vector)
		constexpr inline Vector2D unit () const
		{
			float len = length();
			if (len != 0) {
				return (*this) / len;
			} else {
				return Vector2D(0, 0);
			}
		}
		constexpr inline Vector2D normalized() const
		{
			return unit();
		}
		constexpr inline void normalize()
		{
			*this = unit();
		}

		// Get the orthogonal vector
		constexpr inline Vector2D orthoLeft () const { return Vector2D(-y, x); }
		constexpr inline Vector2D orthoRight () const { return Vector2D(y, -x); }

		// Cross product (the Z component of it)
		constexpr inline T cross (Vector2D param) const { return x * param.y - y * param.x; }

		// Dot product
		constexpr inline T dot (Vector2D param) const { return (x * param.x) + (y * param.y); }

		// Length
		constexpr inline T length () const { return static_cast<T>(std::sqrt(squaredLength())); }
		constexpr inline T len () const { return length(); }
		constexpr inline T manhattanLength() const { return std::abs(x) + std::abs(y); }

		// Squared length, often useful and much faster
		constexpr inline T squaredLength () const { return x*x+y*y; }

		// Projection on another vector
		inline Vector2D projection (Vector2D param) const
		{
			Vector2D unit = param.unit();
			return this->dot(unit) * unit;
		}
		inline T projectionLength (Vector2D param) const
		{
			return this->dot(param.unit());
		}

		// Rounding
		constexpr inline Vector2D floor() const { return Vector2D(std::floor(x), std::floor(y)); }
		constexpr inline Vector2D ceil() const { return Vector2D(std::ceil(x), std::ceil(y)); }
		constexpr inline Vector2D round() const { return Vector2D(std::round(x), std::round(y)); }

		// Gets the angle that this vector is pointing to
		inline U angle () const
		{
			U angle;
			angle.setRadians(atan2(y, x));
			return angle;
		}

		// Rotates vector by an angle
		inline Vector2D rotate (U angle) const
		{
			T cos, sin;
			angle.sincos(sin, cos);
			return Vector2D(x*cos - y*sin, x*sin + y*cos);
		}

		// Rotates vector by sine and cosine
		inline Vector2D rotate (T sine, T cosine) const
		{
			return Vector2D(x*cosine - y*sine, x*sine + y*cosine);
		}

		// Removes the length of the vector along the given axis
		constexpr inline Vector2D neutralize (Vector2D param) const
		{
			return *this - dot(param)*param;
		}

		String toString() const
		{
			return String("(") + x + ", " + y + ")";
		}

		constexpr Vector2D abs() const
		{
			return Vector2D(std::abs(x), std::abs(y));
		}

		constexpr inline static Vector2D<T,U> min(Vector2D<T,U> a, Vector2D<T,U> b)
		{
			return Vector2D<T,U>(std::min(a.x, b.x), std::min(a.y, b.y));
		}

		constexpr inline static Vector2D<T,U> max(Vector2D<T,U> a, Vector2D<T,U> b)
		{
			return Vector2D<T,U>(std::max(a.x, b.x), std::max(a.y, b.y));
		}
	};

	////////////////////
	// Global operators
	template <typename T,class U,typename V>
	inline Vector2D<T,U> operator * (V f, const Vector2D<T,U> &v)
	{
		return Vector2D<T,U>(T(v.x * f),T(v.y * f));
	}

	template <typename T,class U>
	std::ostream& operator<< (std::ostream& ostream, const Vector2D<T, U>& v)
	{
		ostream << "(" << v.x << "," << v.y << ")" ; return ostream;
	}

	////////////
	// Typedefs
	typedef Vector2D<double,Angle<double> > Vector2d;
	typedef Vector2D<> Vector2f;
	typedef Vector2D<int> Vector2i;
	typedef Vector2D<short> Vector2s;
	typedef Vector2D<char> Vector2c;
	typedef Vector2f Position;
	typedef Vector2f Size;
}

namespace std {
	template<>
	struct hash<Halley::Vector2i>
	{
		size_t operator()(const Halley::Vector2i& v) const
		{
			return std::hash<long long>()(*reinterpret_cast<const long long*>(&v));
		}
	};

	template<>
	struct hash<Halley::Vector2f>
	{
		size_t operator()(const Halley::Vector2f& v) const
		{
			return std::hash<long long>()(*reinterpret_cast<const long long*>(&v));
		}
	};
}

max_stars_repo_name: lye/halley

max_issues_repo_name: lye/halley

max_forks_repo_name: lye/halley


text:
/*
  ferrers.c
    Ferrers bar model.
*/

#include <stddef.h>
#include <stdlib.h>
#include <math.h>

#include "utils.h"
#include "potential.h"

#include <gsl/gsl_sf.h>

#include "ferrers.h"

#define Sq(X) ((X) * (X))

double
_lambda(double x_sq, double y_sq, double z_sq, double a_sq, double b_sq, double c_sq)
{
  #define EQS(u) (x_sq / (a_sq + u) + y_sq / (b_sq + u) + z_sq / (c_sq + u) - 1.)
  double J = 0, K = 0., L = (x_sq + y_sq + z_sq);
  for(;EQS(L) > 0;) L *= 2.;//, printf("doubled!\n");

  // solve the eqs using bisection method.
  for(;L - J > 1.e-8;)
    {
      K = (L + J) / 2.;
      //printf("%f, %f, %f\n", J, K, L);
      if(EQS(K) > 0) J = K;
      else L = K;
    }

  return K;

  #undef EQS
}

int
_orbita_potential_ferrers2_evaluate_wijk(double x, double y, double z,
                                         double a, double b, double c,
                                         double * W, int init,
                                         struct orbita_potential_ferrers2_paramset * ws)
{
  double a_sq = a * a, b_sq = b * b, c_sq = c * c,
         x_sq = x * x, y_sq = y * y, z_sq = z * z;

  // outside the bar region?
  if(x_sq / a_sq + y_sq / b_sq + z_sq / c_sq >= 1 || init)
    {
      // using the root-finding eqs. this is NOT the unique positive root.
      // now use bisection method instead.
      /*
      double A = (a_sq + b_sq + c_sq) - (x_sq + y_sq + z_sq),
             B = (a_sq * b_sq + b_sq * c_sq + c_sq * a_sq) -
                 (x_sq * (b_sq + c_sq) + y_sq * (c_sq + a_sq) + z_sq * (a_sq + b_sq)),
             C = (a_sq * b_sq * c_sq) -
                 (x_sq * b_sq * c_sq + a_sq * y_sq * c_sq + a_sq * b_sq * z_sq);

      double A_sq = A * A, B_sq = B * B, C_sq = C * C;
      double A_cb = A_sq * A, B_cb = B_sq * B;
      double CBR_2 = pow(2., 1. / 3.);

      double lambda =
        pow(-2. * A_cb + 3. * sqrt(3.) * sqrt(4. * A_cb * C - A_sq * B_sq - 18. * A * B * C
          + 4. * B_cb + 27. * C_sq) + 9. * A * B - 27. * C, 1. / 3.) / (3. * CBR_2)
          - (CBR_2 * (2187. * B - 729. * A_sq)) / (2187. * pow(-2. * A_cb + 3. * sqrt(3.)
          * sqrt(4. * A_cb * C - A_sq * B_sq - 18. * A * B * C + 4. * B_cb + 27. * C_sq)
          + 9. * A * B - 27. * C, 1. / 3.)) - A / 3.;
      */

      double lambda = _lambda(x_sq, y_sq, z_sq, a_sq, b_sq, c_sq);
      if(init) lambda = 0.;

      double p = asin(sqrt((a_sq - c_sq) / (a_sq + lambda))),
             k = sqrt((a_sq - b_sq) / (a_sq - c_sq));

      /*
      double F = gsl_sf_ellint_F(p, k, GSL_PREC_DOUBLE),
             E = gsl_sf_ellint_E(p, k, GSL_PREC_DOUBLE);
      */

      double F = gsl_sf_ellint_F(p, k, GSL_PREC_SINGLE),
             E = gsl_sf_ellint_E(p, k, GSL_PREC_SINGLE);

      double delta_lambda = sqrt((a_sq + lambda) * (b_sq + lambda) * (c_sq + lambda));

      /* W_000 */ W[ 0] = 2. * F / sqrt(a_sq - c_sq),
      /* W_100 */ W[ 1] = 2. * (F - E) / ((a_sq - b_sq) * sqrt(a_sq - c_sq)),
      /* W_001 */ W[ 3] = (2. / (b_sq - c_sq)) * sqrt((b_sq + lambda) / ((a_sq + lambda) * (c_sq + lambda)))
                          -2. * E / ((b_sq - c_sq) * sqrt(a_sq - c_sq)),
      /* W_010 */ W[ 2] = 2. / delta_lambda - W[1] - W[3],
      /* W_110 */ W[ 4] = (W[2] - W[1]) / (a_sq - b_sq),
      /* W_011 */ W[ 5] = (W[3] - W[2]) / (b_sq - c_sq),
      /* W_101 */ W[ 6] = (W[1] - W[3]) / (c_sq - a_sq),
      /* W_200 */ W[ 7] = (2. / (delta_lambda * (a_sq + lambda)) - W[4] - W[6]) / 3.,
      /* W_020 */ W[ 8] = (2. / (delta_lambda * (b_sq + lambda)) - W[5] - W[4]) / 3.,
      /* W_002 */ W[ 9] = (2. / (delta_lambda * (c_sq + lambda)) - W[6] - W[5]) / 3.,
      /* W_111 */ W[10] = (W[4] - W[5]) / (c_sq - a_sq),
      /* W_120 */ W[11] = (W[8] - W[4]) / (a_sq - b_sq),
      /* W_012 */ W[12] = (W[9] - W[5]) / (b_sq - c_sq),
      /* W_201 */ W[13] = (W[7] - W[6]) / (c_sq - a_sq),
      /* W_210 */ W[14] = (W[7] - W[4]) / (b_sq - a_sq), // NOTE Fixed. Nov. 15, 2015
      /* W_021 */ W[15] = (W[8] - W[5]) / (c_sq - b_sq), //
      /* W_102 */ W[16] = (W[9] - W[6]) / (a_sq - c_sq), //
      /* W_300 */ W[17] = (2. / (delta_lambda * Sq(a_sq + lambda)) - W[14] - W[13]) / 5.,
      /* W_030 */ W[18] = (2. / (delta_lambda * Sq(b_sq + lambda)) - W[15] - W[11]) / 5.,
      /* W_003 */ W[19] = (2. / (delta_lambda * Sq(c_sq + lambda)) - W[16] - W[12]) / 5.;
    }
  else
    {
      int i_ijk;
      for(i_ijk = 0; i_ijk < 20; ++ i_ijk)
        W[i_ijk] = (ws -> W_i)[i_ijk];
    }

  return 0;
}

double
_orbita_potential_ferrers2_psi(const void * par, const double * pos)
{
  struct orbita_potential_ferrers2_paramset * ws =
    (struct orbita_potential_ferrers2_paramset *) par;

  double W[20];
  _orbita_potential_ferrers2_evaluate_wijk(
      pos[0], pos[1], pos[2], ws -> a, ws -> b, ws -> c, W, 0, ws);

  double x_sq = pos[0] * pos[0], y_sq = pos[1] * pos[1], z_sq = pos[2] * pos[2];

  /*
  return (-(ws -> C) / 6.) *
         (
             W[0] - 6. * x_sq * y_sq * z_sq * W[10]
           + x_sq * (x_sq * (3. * W[7] - x_sq * W[17])
           + 3. * (y_sq * (2 * W[4] - y_sq * W[11] - x_sq * W[14]) - W[1]))
           + y_sq * (y_sq * (3. * W[8] - y_sq * W[18])
           + 3. * (z_sq * (2 * W[5] - z_sq * W[12] - y_sq * W[15]) - W[2]))
           + z_sq * (z_sq * (3. * W[9] - z_sq * W[19])
           + 3. * (x_sq * (2 * W[6] - x_sq * W[13] - z_sq * W[16]) - W[3]))
         );
  */

  // This is to test if Del^2 Psi == 4 Pi G rho. Returns rho from W_jkl.
  return ((ws -> a) * (ws -> b) * (ws -> c) * (ws -> rho) / 2.) *
          (
                  (W[1] + W[2] + W[3])
            - 2. * ((y_sq + x_sq) * W[4] + (z_sq + y_sq) * W[5] + (x_sq + z_sq) * W[6])
            - 6. * (x_sq * W[7] + y_sq * W[8] + z_sq * W[9])
            + 2. * (x_sq * y_sq + y_sq * z_sq + z_sq * x_sq) * W[10]
            + (y_sq * y_sq + 6. * x_sq * y_sq) * W[11]
            + (z_sq * z_sq + 6. * y_sq * z_sq) * W[12]
            + (x_sq * x_sq + 6. * z_sq * x_sq) * W[13]
            + (6. * x_sq * y_sq + x_sq * x_sq) * W[14]
            + (6. * y_sq * z_sq + y_sq * y_sq) * W[15]
            + (6. * z_sq * x_sq + z_sq * z_sq) * W[16]
            + 5. * (x_sq * x_sq * W[17] + y_sq * y_sq * W[18] + z_sq * z_sq * W[19])
          );
}

int
_orbita_potential_ferrers2_f(const void * par,
                             const double * pos, double * F)
{
  struct orbita_potential_ferrers2_paramset * ws =
    (struct orbita_potential_ferrers2_paramset *) par;

  double W[20];
  _orbita_potential_ferrers2_evaluate_wijk(
    pos[0], pos[1], pos[2], ws -> a, ws -> b, ws -> c, W, 0, ws);

  double x_sq = pos[0] * pos[0], y_sq = pos[1] * pos[1], z_sq = pos[2] * pos[2]; //

  F[0] = -pos[0] * (ws -> C) * //
          (
            W[1] + x_sq * (x_sq * W[17] + 2. * (y_sq * W[14] - W[7]))
                 + y_sq * (y_sq * W[11] + 2. * (z_sq * W[10] - W[4]))
                 + z_sq * (z_sq * W[16] + 2. * (x_sq * W[13] - W[6]))
          ),
  F[1] = -pos[1] * (ws -> C) * //
          (
            W[2] + x_sq * (x_sq * W[14] + 2. * (y_sq * W[11] - W[4]))
                 + y_sq * (y_sq * W[18] + 2. * (z_sq * W[15] - W[8]))
                 + z_sq * (z_sq * W[12] + 2. * (x_sq * W[10] - W[5]))
          ),
  F[2] = -pos[2] * (ws -> C) * //
          (
            W[3] + x_sq * (x_sq * W[13] + 2. * (y_sq * W[10] - W[6]))
                 + y_sq * (y_sq * W[15] + 2. * (z_sq * W[12] - W[5]))
                 + z_sq * (z_sq * W[19] + 2. * (x_sq * W[16] - W[9]))
          );

  return 0;
}

double
_orbita_potential_ferrers2_rho(const void * par, const double * pos)
{
  struct orbita_potential_ferrers2_paramset * ws =
    (struct orbita_potential_ferrers2_paramset *) par;

  double m_sq = Sq(pos[0]) / Sq(ws -> a) + Sq(pos[1]) / Sq(ws -> b) + Sq(pos[2]) / Sq(ws -> c);
  if(m_sq < 1.) return (ws -> rho) * Sq(1 - m_sq);
  else return 0.;
}

double
_orbita_potential_ferrers2_Mr(const void * par, const double R)
{
  orbita_err("Method not supported.");
  return 0.;
}

double
_orbita_potential_ferrers2_Mt(const void * par)
{
  struct orbita_potential_ferrers2_paramset * ws =
    (struct orbita_potential_ferrers2_paramset *) par;
  return (ws -> rho) * (ws -> a) * (ws -> b) * (ws -> c) * (32. * Pi / 105.);
}

int
_orbita_potential_ferrers2_kill(struct orbita_potential * psi)
{
  free(psi -> param);
  return 0;
}

struct orbita_potential *
orbita_potential_ferrers2(double rho, double a, double b, double c)
{
  struct orbita_potential * psi =
    (struct orbita_potential *) malloc(sizeof(struct orbita_potential));
  struct orbita_potential_ferrers2_paramset * ws =
    (struct orbita_potential_ferrers2_paramset *)
    malloc(sizeof(struct orbita_potential_ferrers2_paramset));
  psi -> param = (void *)ws;

  // initialize w_i
  _orbita_potential_ferrers2_evaluate_wijk(0., 0., 0., a, b, c, ws -> W_i, 1, ws);
  ws -> C = 2. * Pi * rho * a * b * c;
  ws -> a = a, ws -> b = b, ws -> c = c, ws -> rho = rho;

  // set other parameters
  psi -> symmetry     = 0,
  psi -> coordinates  = ORBITA_COORD_CARTESIAN;

  psi -> rho = & _orbita_potential_ferrers2_rho,
  psi -> psi = & _orbita_potential_ferrers2_psi,
  psi -> f   = & _orbita_potential_ferrers2_f,
  psi -> Mr  = & _orbita_potential_ferrers2_Mr,
  psi -> Mt  = & _orbita_potential_ferrers2_Mt;

  psi -> is_variable = 0,
  psi -> init_param  = NULL,
  psi -> param_func  = NULL;

  psi -> kill = & _orbita_potential_ferrers2_kill;

  psi -> is_composite = 0;

  return psi;
}

max_stars_repo_name: shiaki/orbita

max_issues_repo_name: shiaki/orbita

max_forks_repo_name: shiaki/orbita


text:
/*==================BEGIN ASF DOCUMENTATION==================*/
/*
ABOUT EDITING THIS DOCUMENTATION:
If you wish to edit the documentation for this program, you need to change the
following defines.
*/

#define ASF_NAME_STRING \
"asf_export"

#define ASF_USAGE_STRING \
"   "ASF_NAME_STRING" [-format <output_format>] [-byte <sample mapping option>]\n"\
"              [-rgb <red> <green> <blue>] [-band <band_id | all>]\n"\
"              [-lut <look up table file>] [-truecolor] [-falsecolor]\n"\
"              [-log <log_file>] [-quiet] [-license] [-version] [-help]\n"\
"              <in_base_name> <out_full_name>\n"

#define ASF_DESCRIPTION_STRING \
"   This program ingests ASF internal format data and exports said data to a\n"\
"   number of graphics file formats (TIFF/GEOTIFF, JPEG, PGM, PNG, POlSARPRO,\n"\
"   HDF5 and netCDF).\n"\
"   If the input data was geocoded and the ouput format supports geocoding,\n"\
"   that information will be included.  Optionally, you may apply look-up tables,\n"\
"   assign color bands (-rgb, -truecolor, -falsecolor).\n"

#define ASF_INPUT_STRING \
"   A file set in the ASF internal data format.\n"

#define ASF_OUTPUT_STRING \
"   The converted data in the output file, with the requested format.\n"

#define ASF_OPTIONS_STRING \
"   -format <format>\n"\
"        Format to export to. Must be one of the following:\n"\
"            tiff      - Tagged Image File Format, with byte valued pixels\n"\
"            geotiff   - GeoTIFF file, with floating point or byte valued pixels\n"\
"            jpeg      - Lossy compressed image, with byte valued pixels\n"\
"            pgm       - Portable graymap image, with byte valued pixels\n"\
"            png       - Portable network graphic, with byte valued pixels\n"\
"            polsarpro - Flat binary floating point files in PolSARPro format\n"\
"            hdf5      - HDF5 format, with floating point or byte valued pixels\n"\
"            netcdf    - netCDF format compliant to the CF conventions\n"\
"   NOTE: When exporting to a GeoTIFF format file, all map-projection\n"\
"         information is included in GeoKeys as specified in the GeoTIFF\n"\
"         standard.  The other graphics file formats do not support the\n"\
"         storing of map-projection parameters in the output file.  If you\n"\
"         wish to maintain the map-projection and/or georeference (corner\n"\
"         point) information in the output, you should choose the GeoTIFF\n"\
"         output format.\n\n"\
"   NOTE: When exporting to a GeoTIFF format file, the data format (floating\n"\
"         point, byte, 16-bit integer, etc) will be maintained.  Many viewers\n"\
"         cannot view non-integer data.  Leaving the data format the same as\n"\
"         the original produces the most accurate export, but remapping it to\n"\
"         byte range (0-255) with the -byte option will result in the greatest\n"\
"         compatibility with viewers and GIS software packages.\n\n"\
"   -byte <sample mapping option>\n"\
"        Converts output image to byte using the following options:\n"\
"            truncate\n"\
"                values less than 0 are mapped to 0, values greater than 255\n"\
"                are mapped to 255, and values in between are converted to\n"\
"                whole numbers (the fractional part of the values are truncated,\n"\
"                not rounded.)\n"\
"            minmax\n"\
"                determines the minimum and maximum values of the input image\n"\
"                and linearly maps those values to the byte range of 0 to 255.\n"\
"                The remapping is accomplished using real (floating point)\n"\
"                numbers then the result is converted to a whole number by\n"\
"                truncating the fractional part.\n"\
"            sigma\n"\
"                determines the mean and standard deviation of an image and\n"\
"                defines a range of +/- two standard deviations (2-sigma)\n"\
"                around the mean value, and maps this buffer to the byte range\n"\
"                0 to 255 as described for minmax above.  The range limits are\n"\
"                adjusted if the either of the 2-sigma range limits lie outside\n"\
"                the range of the original values.  As with the other remapping\n"\
"                methods, the calculation of values are made with real numbers\n"\
"                and the result is converted to a whole number by truncating\n"\
"                any fractional part.\n"\
"            histogram_equalize\n"\
"                develops a look-up table by integrating the (no-adaptation)\n"\
"                whole-image histogram and then normalizing the result to the\n"\
"                0-255 range.  The result is that areas of low contrast, i.e.\n"\
"                flat topography, have a stronger contrast expansion applied\n"\
"                while areas of high contrast, i.e. non-flat topography, will\n"\
"                receive less contrast expansion.  Histogram equalization\n"\
"                is a useful transform for making hard-to-see detail more\n"\
"                visible for the the viewer but is likewise a nonlinear\n"\
"                transform that results in minor (apparent) topography shifts\n"\
"                within the image.  Since shifts occur the most in areas where\n"\
"                the contrast is expanded the most, i.e. flat topography, and\n"\
"                less in areas of more interesting topography, the pragmatic\n"\
"                conclusion is that the nonlinear shifts are quite\n"\
"                insignificant for the majority of users ...only important\n"\
"                if performing precision geography measurements or overlays.\n"\
"   -rgb <red> <green> <blue>\n"\
"        Converts output image into a color RGB image.\n"\
"        <red>, <green>, and <blue> specify which band (channel)\n"\
"        is assigned to color planes red, green, or blue,\n"\
"        ex) '-rgb HH VH VV', or '-rgb 3 2 1'.  If the word 'ignore' is\n"\
"        provided as one or more bands, then the associated color plane\n"\
"        will not be exported, e.g. '-rgb ignore 2 1' will result in an\n"\
"        RGB file that has a zero in each RGB pixel's red component, band 2\n"\
"        assigned to the green channel, and band 1 assigned to the blue.\n"\
"        The result will be an image with only greens and blues in it.\n"\
"        Currently implemented for GeoTIFF, TIFF, JPEG and PNG.\n"\
"        Cannot be used together with the -band option.\n"\
"   -lut <look up table file>\n"\
"        Applies a color look up table to the image while exporting.\n"\
"        Only allowed for single-band images.  Images must contain byte (8-bit)\n"\
"        data (except for data deriving from PolSARpro classifications.) Some\n"\
"        look-up table files are in the look_up_tables subdirectory in\n"\
"        the asf_tools share directory.  The tool will look in\n"\
"        this directory for the specified file if it isn't found\n"\
"        in the current directory.\n"\
"   -truecolor\n"\
"        For 3 or 4 band optical satellite images where the first band is the\n"\
"        the blue band, the second green, and the third red.  This option will\n"\
"        export the third band as the red element, the second band as the green\n"\
"        element, and the first band as the blue element.  Performs a 2-sigma\n"\
"        constrast expansion on each individual band during the export (similar\n"\
"        to most GIS software packages.)  To export a true-color image WITHOUT\n"\
"        the contrast expansion associated with the truecolor option, use the\n"\
"        -rgb option instead.  The rgb option will directly assign the\n"\
"        available bands, unaltered, to the RGB color channels in the output\n"\
"        file, e.g.\n\n"\
"          'asf_export -rgb 03 02 01 <infile> <outfile>'.\n\n"\
"        Only allowed for multi-band images with 3 or more bands.\n"\
"        The truecolor option cannot be used together with any of the\n"\
"        following options: -rgb, -band, or -falsecolor.\n"\
"   -falsecolor\n"\
"        For 4 band optical satellite images where the second band is green,\n"\
"        the third red, and the fourth is the near-infrared band.  Exports\n"\
"        the fourth (IR) band as the red element, the third band as the green\n"\
"        element, and the second band as the blue element.  Performs a 2-sigma\n"\
"        constrast expansion on the individual bands during the export.  To\n"\
"        export a falsecolor image WITHOUT the contrast expansion, use the\n"\
"        -rgb flag to directly assign the available bands to the RGB color\n"\
"        channels, e.g.\n\n"\
"          'asf_export -rgb 04 03 02 <infile> <outfile>'.\n\n"\
"        Only allowed for multi-band images with 4 bands.\n"\
"        Cannot be used together with any of the following: -rgb, -band,\n"\
"        or -truecolor.\n"\
"   -band <band_id | all>\n"\
"        If the data contains multiple data files, one for each band (channel)\n"\
"        then export the band identified by 'band_id' (only).  If 'all' is\n"\
"        specified rather than a band_id, then export all available bands into\n"\
"        individual files, one for each band.  Default is '-band all'.\n"\
"        Cannot be chosen together with the -rgb option.\n"\
"   -log <logFile>\n"\
"        Output will be written to a specified log file.\n"\
"   -quiet\n"\
"        Supresses all non-essential output.\n"\
"   -license\n"\
"        Print copyright and license for this software then exit.\n"\
"   -version\n"\
"        Print version and copyright then exit.\n"\
"   -help\n"\
"        Print a help page and exit.\n"

#define ASF_EXAMPLES_STRING \
"   To export to the default GeoTIFF format from file1.img and file1.meta\n"\
"   to file1.tif:\n"\
"        example> "ASF_NAME_STRING" file1 file1\n\n"\
"   NOTE: When exporting to a GeoTIFF format file, all map-projection\n"\
"         information is included in GeoKeys as specified in the GeoTIFF\n"\
"         standard.\n\n"\
"   NOTE: When exporting to a GeoTIFF format file, the data format (floating\n"\
"         point, byte, 16-bit integer, etc) will be maintained.  Many viewers\n"\
"         cannot view non-integer data.  Leaving the data format the same as\n"\
"         the original produces the most accurate export, but remapping it to\n"\
"         byte range (0-255) with the -byte option will result in the greatest\n"\
"         compatibility with viewers and GIS software packages.\n"\
"\n"\
"   To export to file2.jpg in the jpeg format:\n"\
"        example> "ASF_NAME_STRING" -format jpeg file1 file2\n"\
"\n"

#define ASF_LIMITATIONS_STRING \
"   Currently supports ingest of ASF format floating point and byte data (only).\n"\
"\n"\
"   Floating-point image formats (i.e., geotiff) are not supported in many\n"\
"   image viewing programs.\n"

#define ASF_SEE_ALSO_STRING \
"   asf_mapready, asf_import\n"

/*===================END ASF DOCUMENTATION===================*/

#include <ctype.h>
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <unistd.h>
#include <limits.h>

#include <cla.h>
#include <envi.h>
#include <esri.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_statistics.h>

#include <asf.h>
#include <asf_endian.h>
#include <asf_meta.h>
#include <asf_raster.h>
#include <asf_export.h>
#include <asf_contact.h>
#include <asf_license.h>

// Local prototypes
int is_numeric(char *str);

// Print minimalistic usage info & exit
static void print_usage(void)
{
  asfPrintStatus("\n"
      "Usage:\n"
      ASF_USAGE_STRING
      "\n");
  exit(EXIT_FAILURE);
}

// Print the help info & exit
static void print_help(void)
{
  asfPrintStatus(
      "\n"
      "Tool name:\n   " ASF_NAME_STRING "\n\n"
      "Usage:\n" ASF_USAGE_STRING "\n"
      "Description:\n" ASF_DESCRIPTION_STRING "\n"
      "Input:\n" ASF_INPUT_STRING "\n"
      "Output:\n"ASF_OUTPUT_STRING "\n"
      "Options:\n" ASF_OPTIONS_STRING "\n"
      "Examples:\n" ASF_EXAMPLES_STRING "\n"
      "Limitations:\n" ASF_LIMITATIONS_STRING "\n"
      "See also:\n" ASF_SEE_ALSO_STRING "\n"
      "Contact:\n" ASF_CONTACT_STRING "\n"
      "Version:\n   " SVN_REV " (part of " TOOL_SUITE_NAME " " MAPREADY_VERSION_STRING ")\n\n");
  exit(EXIT_FAILURE);
}

int
checkForOption (char *key, int argc, char *argv[])
{
  int ii = 0;
  while ( ii < argc ) {
    if ( strmatch (key, argv[ii]) )
      return ii;
    ++ii;
  }
  return FLAG_NOT_SET;
}

static const char *sigma_str(int with_sigma) {
    return with_sigma ? "w/sigma" : "";
}

// Main program body.
int
main (int argc, char *argv[])
{
  output_format_t format = 0;
  meta_parameters *md;
  char *in_base_name, *output_name;
  char **band_names=NULL;
  int rgb=0;
  int true_color;
  int false_color;
  int num_bands_found;
  int ignored[3] = {0, 0, 0};
  int num_ignored = 0;

  in_base_name = (char *) MALLOC(sizeof(char)*255);
  output_name = (char *) MALLOC(sizeof(char)*255);

/**********************BEGIN COMMAND LINE PARSING STUFF**********************/
  // Command line input goes in it's own structure.
  command_line_parameters_t command_line;
  strcpy (command_line.format, "");
  command_line.size = NO_MAXIMUM_OUTPUT_SIZE;
  strcpy (command_line.in_data_name, "");
  strcpy (command_line.in_meta_name, "");
  strcpy (command_line.output_name, "");
  command_line.verbose = FALSE;
  command_line.quiet = FALSE;
  strcpy (command_line.leader_name, "");
  strcpy (command_line.cal_params_file, "");
  strcpy (command_line.cal_comment, "");
  command_line.sample_mapping = 0;
  strcpy(command_line.red_channel, "");
  strcpy(command_line.green_channel, "");
  strcpy(command_line.blue_channel, "");
  strcpy(command_line.band, "");
  strcpy(command_line.look_up_table_name, "");

  int formatFlag, logFlag, quietFlag, byteFlag, rgbFlag, bandFlag, lutFlag;
  int truecolorFlag, falsecolorFlag;
  int needed_args = 3;  //command & argument & argument
  int ii;
  char sample_mapping_string[25];

  //Check to see which options were specified
  if (   (checkForOption("--help", argc, argv) != FLAG_NOT_SET)
      || (checkForOption("-h", argc, argv) != FLAG_NOT_SET)
      || (checkForOption("-help", argc, argv) != FLAG_NOT_SET) ) {
      print_help();
  }
  get_asf_share_dir_with_argv0(argv[0]);
  handle_license_and_version_args(argc, argv, ASF_NAME_STRING);

  formatFlag = checkForOption ("-format", argc, argv);
  logFlag = checkForOption ("-log", argc, argv);
  quietFlag = checkForOption ("-quiet", argc, argv);
  byteFlag = checkForOption ("-byte", argc, argv);
  rgbFlag = checkForOption ("-rgb", argc, argv);
  bandFlag = checkForOption ("-band", argc, argv);
  lutFlag = checkForOption ("-lut", argc, argv);
  truecolorFlag = checkForOption("-truecolor", argc, argv);
  falsecolorFlag = checkForOption("-falsecolor", argc, argv);

  if ( formatFlag != FLAG_NOT_SET ) {
    needed_args += 2;           // Option & parameter.
  }
  if ( quietFlag != FLAG_NOT_SET ) {
    needed_args += 1;           // Option & parameter.
  }
  if ( logFlag != FLAG_NOT_SET ) {
    needed_args += 2;           // Option & parameter.
  }
  if ( byteFlag != FLAG_NOT_SET ) {
    needed_args += 2;           // Option & parameter.
  }
  if ( rgbFlag != FLAG_NOT_SET ) {
    needed_args += 4;           // Option & 3 parameters.
  }
  if ( bandFlag != FLAG_NOT_SET ) {
    needed_args += 2;           // Option & parameter.
  }
  if ( lutFlag != FLAG_NOT_SET ) {
    needed_args += 2;           // Option & parameter.
  }
  if ( truecolorFlag != FLAG_NOT_SET ) {
    needed_args += 1;           // Option only
  }
  if ( falsecolorFlag != FLAG_NOT_SET ) {
    needed_args += 1;           // Option only
  }

  if ( argc != needed_args ) {
    print_usage ();                   // This exits with a failure.
  }

  // We also need to make sure the last three options are close to
  // what we expect.
  if ( argv[argc - 1][0] == '-' || argv[argc - 2][0] == '-' ) {
    print_usage (); // This exits with a failure.
  }

  // Make sure any options that have parameters are followed by
  // parameters (and not other options) Also make sure options'
  // parameters don't bleed into required arguments.
  if ( formatFlag != FLAG_NOT_SET ) {
    if ( argv[formatFlag + 1][0] == '-' || formatFlag >= argc - 3 ) {
      print_usage ();
    }
  }
  if ( byteFlag != FLAG_NOT_SET ) {
    if ( argv[byteFlag + 1][0] == '-' || byteFlag >= argc - 3 ) {
      print_usage ();
    }
  }
  if ( rgbFlag != FLAG_NOT_SET ) {
    if (( argv[rgbFlag + 1][0] == '-' && argv[rgbFlag + 2][0] == '-' &&
          argv[rgbFlag + 3][0] == '-' ) || rgbFlag >= argc - 5 ) {
      print_usage ();
    }
  }
  if ( bandFlag != FLAG_NOT_SET ) {
    if ( argv[bandFlag + 1][0] == '-' || bandFlag >= argc - 3 ) {
      print_usage ();
    }
  }
  if ( lutFlag != FLAG_NOT_SET ) {
    if ( argv[lutFlag + 1][0] == '-' || lutFlag >= argc - 3 ) {
      print_usage ();
    }
  }
  if ( logFlag != FLAG_NOT_SET ) {
    if ( argv[logFlag + 1][0] == '-' || logFlag >= argc - 3 ) {
      print_usage ();
    }
  }

  // Make sure there are no flag incompatibilities
  if ( (rgbFlag != FLAG_NOT_SET           &&
        (bandFlag != FLAG_NOT_SET         ||
         truecolorFlag != FLAG_NOT_SET    ||
         falsecolorFlag != FLAG_NOT_SET))
     ||
       (bandFlag != FLAG_NOT_SET         &&
        (rgbFlag != FLAG_NOT_SET         ||
         truecolorFlag != FLAG_NOT_SET   ||
         falsecolorFlag != FLAG_NOT_SET))
     ||
       (truecolorFlag != FLAG_NOT_SET    &&
        (bandFlag != FLAG_NOT_SET        ||
         rgbFlag != FLAG_NOT_SET         ||
         falsecolorFlag != FLAG_NOT_SET))
     ||
       (falsecolorFlag != FLAG_NOT_SET   &&
        (bandFlag != FLAG_NOT_SET        ||
         truecolorFlag != FLAG_NOT_SET   ||
         rgbFlag != FLAG_NOT_SET))
     )
  {
    asfPrintWarning("The following options may only be used one at a time:\n"
        "    %s\n    %s\n    %s\n    %s\n    %s\n    %s\n",
        "-rgb", "-truecolor", "-falsecolor", "-band");
    print_help();
  }
  if ( (rgbFlag != FLAG_NOT_SET         ||
        truecolorFlag != FLAG_NOT_SET   ||
        falsecolorFlag != FLAG_NOT_SET) &&
        lutFlag != FLAG_NOT_SET
     )
    asfPrintError("Look up table option can only be used on single-band "
          "images.\n");

  if( logFlag != FLAG_NOT_SET ) {
    strcpy(logFile, argv[logFlag+1]);
  }
  else {
    sprintf(logFile, "tmp%i.log", (int)getpid());
  }
  logflag = TRUE; // Since we always log, set the old school logflag to true
  fLog = FOPEN (logFile, "a");

  // Set old school quiet flag (for use in our libraries)
  quietflag = ( quietFlag != FLAG_NOT_SET ) ? TRUE : FALSE;

  // We're good enough at this point... print the splash screen.
  asfSplashScreen (argc, argv);

  // Grab the input and output name
  strcpy (in_base_name, argv[argc - 2]);
  strcpy (output_name, argv[argc - 1]);
  strcpy (command_line.output_name, output_name);

  // If user added ".img", strip it.
  char *ext = findExt(in_base_name);
  if (ext && strcmp(ext, ".img") == 0) *ext = '\0';

  // Set default output type
  if( formatFlag != FLAG_NOT_SET ) {
    strcpy (command_line.format, argv[formatFlag + 1]);
  }
  else {
    // Default behavior: produce a geotiff.
    strcpy (command_line.format, "geotiff");
  }

  // Compose input metadata name
  strcpy (command_line.in_meta_name, in_base_name);
  strcat (command_line.in_meta_name, ".meta");

  // for some validation, need the metadata
  md = meta_read (command_line.in_meta_name);

  // Convert the string to upper case.
  for ( ii = 0 ; ii < strlen (command_line.format) ; ++ii ) {
    command_line.format[ii] = toupper (command_line.format[ii]);
  }
  if (strcmp (command_line.format, "PGM") == 0 &&
      (rgbFlag != FLAG_NOT_SET      ||
      truecolorFlag != FLAG_NOT_SET ||
      falsecolorFlag != FLAG_NOT_SET)
     )
  {
    asfPrintWarning("Greyscale PGM output is not compatible with color options:\n"
                    "(RGB, True Color, False Color, color look-up tables, etc\n)"
                    "...Defaulting to producing separate greyscale PGM files for available band.\n");
    rgbFlag = FLAG_NOT_SET;
    truecolorFlag = FLAG_NOT_SET;
    falsecolorFlag = FLAG_NOT_SET;
  }

  // Set the default byte scaling mechanisms
  if (md->optical) {
     // for optical data, default sample mapping is NONE
      command_line.sample_mapping = NONE;
  }
  // for other data, default is based on the output type
  else if (strcmp (command_line.format, "TIFF") == 0 ||
           strcmp (command_line.format, "TIF")  == 0 ||
           strcmp (command_line.format, "JPEG") == 0 ||
           strcmp (command_line.format, "JPG")  == 0 ||
           strcmp (command_line.format, "PNG")  == 0 ||
           strcmp (command_line.format, "PGM")  == 0 ||
	   strcmp (command_line.format, "PNG_ALPHA") == 0 ||
	   strcmp (command_line.format, "PNG_GE") == 0)
  {
    command_line.sample_mapping = SIGMA;
  }
  else if (strcmp (command_line.format, "GEOTIFF") == 0) {
    command_line.sample_mapping = NONE;
  }

  if ( quietFlag != FLAG_NOT_SET )
    command_line.quiet = TRUE;
  else
    command_line.quiet = FALSE;

  // Set rgb combination
  if ( rgbFlag != FLAG_NOT_SET ) {
    int i;

    for (i=0, num_ignored = 0; i<3; i++) {
      ignored[i] = strncmp("IGNORE", uc(argv[rgbFlag + i + 1]), 6) == 0 ? 1 : 0;
      num_ignored += ignored[i] ? 1 : 0;
    }
    asfRequire(num_ignored < 3,
               "Cannot ignore all bands.  Exported image would be blank.\n");

    strcpy (command_line.red_channel, ignored[0] ? "Ignored" : argv[rgbFlag + 1]);
    strcpy (command_line.green_channel, ignored[1] ? "Ignored" : argv[rgbFlag + 2]);
    strcpy (command_line.blue_channel, ignored[2] ? "Ignored" : argv[rgbFlag + 3]);

    // Check to see if the bands are numeric and in range
    int r_channel = atoi(command_line.red_channel);
    int g_channel = atoi(command_line.green_channel);
    int b_channel = atoi(command_line.blue_channel);

    /////////// Numeric channel case ////////////
    // Remove trailing non-numeric characters from the channel number
    // string and pad front end nicely with a zero
    if (!ignored[0] && is_numeric(command_line.red_channel) &&
        r_channel >= 1 && r_channel <= MAX_BANDS) {
      sprintf(command_line.red_channel, "%02d", atoi(command_line.red_channel));
    }
    if (!ignored[1] && is_numeric(command_line.green_channel) &&
        g_channel >= 1 && g_channel <= MAX_BANDS) {
      sprintf(command_line.green_channel, "%02d", atoi(command_line.green_channel));
    }
    if (!ignored[2] && is_numeric(command_line.blue_channel) &&
        b_channel >= 1 && b_channel <= MAX_BANDS) {
      sprintf(command_line.blue_channel, "%02d", atoi(command_line.blue_channel));
    }
  }

  // Set up the bands for true or false color optical data
  true_color = false_color = 0;
  int with_sigma = FALSE;
  if (truecolorFlag != FLAG_NOT_SET || falsecolorFlag != FLAG_NOT_SET) {
    int ALOS_optical = (md->optical && strncmp(md->general->sensor, "ALOS", 4) == 0) ? 1 : 0;
    if (md->optical && truecolorFlag != FLAG_NOT_SET) {
      if (ALOS_optical) {
        with_sigma = TRUE;
        strcpy(command_line.red_channel,   "03");
        strcpy(command_line.green_channel, "02");
        strcpy(command_line.blue_channel,  "01");
        true_color = 1;
        asfPrintStatus("Applying True Color contrast expansion to following channels:");
      }
      else {
        char **bands = extract_band_names(md->general->bands, 3);
        asfRequire(bands != NULL,
                   "-truecolor option specified for non-true color optical image.\n");

        asfPrintWarning("Attempting to use the -truecolor option with non-ALOS\n"
            "optical data.\n");
        strcpy(command_line.red_channel, bands[2]);
        strcpy(command_line.green_channel, bands[1]);
        strcpy(command_line.blue_channel, bands[0]);
        int i;
        for (i=0; i<3; i++) {
          FREE(bands[i]);
        }
        FREE(bands);
      }
    }
    if (md->optical && falsecolorFlag != FLAG_NOT_SET) {
      if (ALOS_optical) {
        with_sigma = TRUE;
        strcpy(command_line.red_channel,   "04");
        strcpy(command_line.green_channel, "03");
        strcpy(command_line.blue_channel,  "02");
        false_color = 1;
        asfPrintStatus("Applying False Color contrast expansion to the following channels:");
      }
      else {
        char **bands = extract_band_names(md->general->bands, 4);
        asfRequire(bands != NULL,
                   "-falsecolor option specified for an optical image with fewer than 4 bands.\n");

        asfPrintWarning("Attempting to use the -falsecolor option with non-ALOS\n"
            "optical data.\n");
        strcpy(command_line.red_channel, bands[3]);
        strcpy(command_line.green_channel, bands[2]);
        strcpy(command_line.blue_channel, bands[1]);
        int i;
        for (i=0; i<3; i++) {
          FREE(bands[i]);
        }
        FREE(bands);
      }
    }
    if (!ALOS_optical && !md->optical) {
      asfPrintError("-truecolor or -falsecolor option selected with non-optical data\n");
    }
  }

  if (rgbFlag != FLAG_NOT_SET ||
     truecolorFlag != FLAG_NOT_SET ||
     falsecolorFlag != FLAG_NOT_SET)
  {
    char red_band[16], green_band[16], blue_band[16];

    asfPrintStatus("\nRed channel  : %s %s\n", command_line.red_channel, sigma_str(with_sigma));
    asfPrintStatus("Green channel: %s %s\n", command_line.green_channel, sigma_str(with_sigma));
    asfPrintStatus("Blue channel : %s %s\n\n", command_line.blue_channel, sigma_str(with_sigma));

    if (is_numeric(command_line.red_channel) &&
        is_numeric(command_line.green_channel) &&
        is_numeric(command_line.blue_channel))
    {
      sprintf(red_band, "%02d", atoi(command_line.red_channel));
      sprintf(green_band, "%02d", atoi(command_line.green_channel));
      sprintf(blue_band, "%02d", atoi(command_line.blue_channel));
      band_names = find_bands(in_base_name, rgbFlag,
                              red_band,
                              green_band,
                              blue_band,
                              &num_bands_found);
    }
    else {
      band_names = find_bands(in_base_name, rgbFlag,
                              command_line.red_channel,
                              command_line.green_channel,
                              command_line.blue_channel,
                              &num_bands_found);
    }
  }

  // Set band
  if ( bandFlag != FLAG_NOT_SET) {
    strcpy (command_line.band, argv[bandFlag + 1]);
    band_names = find_single_band(in_base_name, command_line.band,
                  &num_bands_found);
  }
  else if (rgbFlag == FLAG_NOT_SET &&
          truecolorFlag == FLAG_NOT_SET &&
          falsecolorFlag == FLAG_NOT_SET &&
          bandFlag == FLAG_NOT_SET) {
    bandFlag=1; // For proper messaging to the user
    strcpy (command_line.band, "all");
    band_names = find_single_band(in_base_name, command_line.band,
                                  &num_bands_found);
  }

  // Read look up table name
  if ( lutFlag != FLAG_NOT_SET) {
    strcpy(command_line.look_up_table_name, argv[lutFlag + 1]);
    rgb = 1;
  }

  // Set scaling mechanism
  if ( byteFlag != FLAG_NOT_SET ) {
    strcpy (sample_mapping_string, argv[byteFlag + 1]);
    for ( ii = 0; ii < strlen(sample_mapping_string); ii++) {
      sample_mapping_string[ii] = toupper (sample_mapping_string[ii]);
    }
    if ( strcmp (sample_mapping_string, "TRUNCATE") == 0 )
      command_line.sample_mapping = TRUNCATE;
    else if ( strcmp(sample_mapping_string, "MINMAX") == 0 )
      command_line.sample_mapping = MINMAX;
    else if ( strcmp(sample_mapping_string, "SIGMA") == 0 )
      command_line.sample_mapping = SIGMA;
    else if ( strcmp(sample_mapping_string, "HISTOGRAM_EQUALIZE") == 0 )
      command_line.sample_mapping = HISTOGRAM_EQUALIZE;
    else if ( strcmp(sample_mapping_string, "NONE") == 0 ) {
        asfPrintWarning("Sample remapping method (-byte option) is set to NONE\n"
                "which doesn't make sense.  Defaulting to TRUNCATE...\n");
        command_line.sample_mapping = TRUNCATE;
    }
    else
      asfPrintError("Unrecognized byte scaling method '%s'.\n",
                    sample_mapping_string);
  }

  int is_polsarpro = (md->general->bands && strstr(md->general->bands, "POLSARPRO") != NULL) ? 1 : 0;
  if ( !is_polsarpro               &&
       lutFlag != FLAG_NOT_SET     &&
       bandFlag == FLAG_NOT_SET    &&
       md->general->band_count > 1)
  {
    asfPrintError("Look up tables can only be applied to single band"
          " images\n");
  }

  if ( !is_polsarpro                       &&
       lutFlag != FLAG_NOT_SET             &&
       command_line.sample_mapping == NONE &&
       md->general->data_type != BYTE      &&
       md->general->band_count == 1)
  {
    asfPrintError("Look up tables can only be applied to byte output"
          " images\n");
  }

  // Report what is going to happen
  if (rgbFlag != FLAG_NOT_SET ||
     truecolorFlag != FLAG_NOT_SET ||
     falsecolorFlag != FLAG_NOT_SET)
  {
    if (num_bands_found >= 3) {
      asfPrintStatus("Exporting multiband image ...\n\n");
      rgb = 1;
    }
    else {
      asfPrintError("Not all RGB channels found.\n");
    }
  }
  else if (bandFlag != FLAG_NOT_SET) {
    if (strcmp_case(command_line.band, "ALL") == 0) {
      if (multiband(command_line.format,
		    extract_band_names(md->general->bands, md->general->band_count),
		    md->general->band_count))
	asfPrintStatus("Exporting multiband image ...\n\n");
      else  if (num_bands_found > 1)
        asfPrintStatus("Exporting each band into individual greyscale files ...\n\n");
    }
    else if (num_bands_found == 1) {
      if (lutFlag != FLAG_NOT_SET)
    asfPrintStatus("Exporting band '%s' applying look up table ...\n\n",
               command_line.band);
      else
    asfPrintStatus("Exporting band '%s' as greyscale ...\n\n",
               command_line.band);
    }
    else
      asfPrintError("Band could not be found in the image.\n");
  }
  else if (lutFlag != FLAG_NOT_SET)
    asfPrintStatus("Exporting applying look up table.\n\n");
  else
    asfPrintStatus("Exporting as greyscale.\n\n");

  //If user added ".img", strip it.
  ext = findExt(in_base_name);
  if (ext && strcmp(ext, ".img") == 0) *ext = '\0';
  meta_free(md);

/***********************END COMMAND LINE PARSING STUFF***********************/

  if ( strcmp (command_line.format, "ENVI") == 0 ) {
    format = ENVI;
  }
  else if ( strcmp (command_line.format, "ESRI") == 0 ) {
    format = ESRI;
  }
  else if ( strcmp (command_line.format, "GEOTIFF") == 0 ||
            strcmp (command_line.format, "GEOTIF") == 0) {
    format = GEOTIFF;
  }
  else if ( strcmp (command_line.format, "TIFF") == 0 ||
            strcmp (command_line.format, "TIF") == 0) {
    format = TIF;
  }
  else if ( strcmp (command_line.format, "JPEG") == 0 ||
            strcmp (command_line.format, "JPG") == 0) {
    format = JPEG;
  }
  else if ( strcmp (command_line.format, "PGM") == 0 ) {
    format = PGM;
  }
  else if ( strcmp (command_line.format, "PNG") == 0 ) {
    format = PNG;
  }
  else if ( strcmp (command_line.format, "PNG_ALPHA") == 0 ) {
    format = PNG_ALPHA;
  }
  else if ( strcmp (command_line.format, "PNG_GE") == 0 ) {
    format = PNG_GE;
  }
  else if ( strcmp (command_line.format, "KML") == 0 ) {
    format = KML;
  }
  else if ( strcmp (command_line.format, "POLSARPRO") == 0 ) {
    format = POLSARPRO_HDR;
  }
  else if ( strcmp (command_line.format, "HDF5") == 0 ) {
    format = HDF;
  }
  else if ( strcmp (command_line.format, "NETCDF") == 0 ) {
    format = NC;
  }
  else {
    asfPrintError("Unrecognized output format specified\n");
  }

  /* Complex data generally can't be output into meaningful images, so
     we refuse to deal with it.  */
  /*
  md = meta_read (command_line.in_meta_name);
  asfRequire (   md->general->data_type == BYTE
              || md->general->data_type == INTEGER16
              || md->general->data_type == INTEGER32
              || md->general->data_type == REAL32
              || md->general->data_type == REAL64,
              "Cannot cope with complex data, exiting...\n");

  meta_free (md);
  */

  // Do that exporting magic!
  asf_export_bands(format, command_line.sample_mapping, rgb,
                   true_color, false_color,
                   command_line.look_up_table_name,
                   in_base_name, command_line.output_name, band_names,
                   NULL, NULL);

  // If the user didn't ask for a log file then nuke the one that's been kept
  // since everything has finished successfully
  if (logFlag == FLAG_NOT_SET) {
      fclose (fLog);
      remove(logFile);
  }

  for (ii = 0; ii<num_bands_found; ii++) {
    FREE(band_names[ii]);
  }
  FREE(band_names);
  FREE(in_base_name);
  FREE(output_name);
  exit (EXIT_SUCCESS);
}

int is_numeric(char *str)
{
  char *s = str;
  int numericness = 1;

  while (*s != '\0') {
    if (!isdigit(*s)) numericness = 0;
    s++;
  }

  return numericness;
}


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */

#include "funcs.h"
#include "initial.h"

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_ieee_utils.h>

int main ( int argc, char* argv[] )
{
	gsl_ieee_env_setup () ;			/* read GSL_IEEE_MODE */

	printf("Temperature: %g\n\n", T*HBAR*Delta/BOLTZ) ;

	double beta = 1.0/T ;           	/* Boltzmann factor: beta */
	double omega_1 = gsl_hypot(OMEGA,D) ;   /* omega' */

	struct f_params params;
	params.omega_c = omega_c ;
	params.beta = beta ;
	params.Omega = OMEGA ;
	params.omega_1 = omega_1 ;
	params.alpha = alpha ;

	int status1 = save_integrals ( &params ) ;
	int status2 = save_matrices ( &params ) ;

	/* read the Redfield matrix from a file */
	gsl_matrix* red_m = gsl_matrix_calloc ( 4, 4 ) ;
	int status3 = mat_read ( red_m, "REDFIELD_MATRIX" ) ;

	/* read the CP matrix from a file */
	gsl_matrix* cp_m = gsl_matrix_calloc ( 4, 4 ) ;
	int status4 = mat_read ( cp_m, "CP_MATRIX" ) ;

	/* Hamiltonian generator */
	const double om[] = { 0 , 0 , 0 , omega_1/2 } ;
	gsl_matrix* H = gsl_matrix_calloc ( 4, 4 ) ;
	int status5 = ham_gen ( H, om ) ;

	/*
	 * Find the stationary state by solving the linear systems
	 * and the associated stationary currents
	 */
	gsl_vector* req_red = gsl_vector_calloc (4) ;
	gsl_vector* req_cp = gsl_vector_calloc (4) ;

	printf("REDFIELD DYNAMICS\n") ;
	int status6 = stationary ( red_m , req_red ) ;
	printf("Stationary state: ( %.1f , %.9f , %.9f , %.9f )\n",
			VECTOR(req_red,0), VECTOR(req_red,1),
			VECTOR(req_red,2), VECTOR(req_red,3) ) ;
	printf("Stationary normalized (I/I0) DC current: %.9f\n\n",
			-VECTOR(req_red,3)*OMEGA/omega_1) ;

	printf("CP DYNAMICS\n") ;
	int status7 = stationary ( cp_m , req_cp ) ;
	printf("Stationary state: ( %.1f , %.9f , %.9f , %.9f )\n",
			VECTOR(req_cp,0), VECTOR(req_cp,1),
			VECTOR(req_cp,2), VECTOR(req_cp,3) ) ;
	printf("Stationary normalized (I/I0) DC current: %.9f\n\n",
			-VECTOR(req_cp,3)*OMEGA/omega_1) ;

	/* Save the stationary states into files */
	FILE* f = fopen ( "RED_STATIONARY.dat", "w" ) ;
	gsl_vector_fprintf ( f, req_red, "%.9f" ) ;
	FILE* g = fopen ( "CP_STATIONARY.dat", "w" ) ;
	gsl_vector_fprintf ( g, req_cp, "%.9f" ) ;

	fclose (f) ; fclose (g) ;

	/* polarization */
	double D30 = gsl_matrix_get(cp_m,3,0) ;
	double D33 = gsl_matrix_get(cp_m,3,3) ;
	double pol = -D30/D33 ;
	printf("n-Polarization of the CP dynamics -D30/D33: %.9f\n", pol ) ;

	/* polarization through the formula (for checking) */
	double P = polarization ( &params, omega_c ) ;
	printf("\n") ;
	printf("n-Polarization through the given formula: %.9f\n", P ) ;
	printf("Stationary normalized (I/I0) DC current: %.9f\n", -P*OMEGA/omega_1 ) ;

	/* free memory for matrices */
	gsl_matrix_free(red_m) ;
 	gsl_matrix_free(cp_m) ;

	return status1 + status2 + status3 + status4 + status5 + status6 + status7 ;
}


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:43:55 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -real2hc 12 */

/*
 * This function contains 38 FP additions, 8 FP multiplications,
 * (or, 34 additions, 4 multiplications, 4 fused multiply/add),
 * 18 stack variables, and 24 memory accesses
 */
static const fftw_real K866025403 = FFTW_KONST(+0.866025403784438646763723170752936183471402627);
static const fftw_real K500000000 = FFTW_KONST(+0.500000000000000000000000000000000000000000000);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_real2hc_12(const fftw_real *input, fftw_real *real_output, fftw_real *imag_output, int istride, int real_ostride, int imag_ostride)
{
     fftw_real tmp5;
     fftw_real tmp25;
     fftw_real tmp11;
     fftw_real tmp23;
     fftw_real tmp30;
     fftw_real tmp35;
     fftw_real tmp10;
     fftw_real tmp26;
     fftw_real tmp12;
     fftw_real tmp18;
     fftw_real tmp29;
     fftw_real tmp34;
     fftw_real tmp31;
     fftw_real tmp32;
     ASSERT_ALIGNED_DOUBLE;
     {
	  fftw_real tmp1;
	  fftw_real tmp2;
	  fftw_real tmp3;
	  fftw_real tmp4;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp1 = input[0];
	  tmp2 = input[4 * istride];
	  tmp3 = input[8 * istride];
	  tmp4 = tmp2 + tmp3;
	  tmp5 = tmp1 + tmp4;
	  tmp25 = tmp1 - (K500000000 * tmp4);
	  tmp11 = tmp3 - tmp2;
     }
     {
	  fftw_real tmp19;
	  fftw_real tmp20;
	  fftw_real tmp21;
	  fftw_real tmp22;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp19 = input[9 * istride];
	  tmp20 = input[istride];
	  tmp21 = input[5 * istride];
	  tmp22 = tmp20 + tmp21;
	  tmp23 = tmp19 - (K500000000 * tmp22);
	  tmp30 = tmp19 + tmp22;
	  tmp35 = tmp21 - tmp20;
     }
     {
	  fftw_real tmp6;
	  fftw_real tmp7;
	  fftw_real tmp8;
	  fftw_real tmp9;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp6 = input[6 * istride];
	  tmp7 = input[10 * istride];
	  tmp8 = input[2 * istride];
	  tmp9 = tmp7 + tmp8;
	  tmp10 = tmp6 + tmp9;
	  tmp26 = tmp6 - (K500000000 * tmp9);
	  tmp12 = tmp8 - tmp7;
     }
     {
	  fftw_real tmp14;
	  fftw_real tmp15;
	  fftw_real tmp16;
	  fftw_real tmp17;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp14 = input[3 * istride];
	  tmp15 = input[7 * istride];
	  tmp16 = input[11 * istride];
	  tmp17 = tmp15 + tmp16;
	  tmp18 = tmp14 - (K500000000 * tmp17);
	  tmp29 = tmp14 + tmp17;
	  tmp34 = tmp16 - tmp15;
     }
     real_output[3 * real_ostride] = tmp5 - tmp10;
     imag_output[3 * imag_ostride] = tmp29 - tmp30;
     tmp31 = tmp5 + tmp10;
     tmp32 = tmp29 + tmp30;
     real_output[6 * real_ostride] = tmp31 - tmp32;
     real_output[0] = tmp31 + tmp32;
     {
	  fftw_real tmp37;
	  fftw_real tmp38;
	  fftw_real tmp33;
	  fftw_real tmp36;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp37 = tmp34 + tmp35;
	  tmp38 = tmp11 + tmp12;
	  imag_output[2 * imag_ostride] = K866025403 * (tmp37 - tmp38);
	  imag_output[4 * imag_ostride] = K866025403 * (tmp38 + tmp37);
	  tmp33 = tmp25 - tmp26;
	  tmp36 = K866025403 * (tmp34 - tmp35);
	  real_output[5 * real_ostride] = tmp33 - tmp36;
	  real_output[real_ostride] = tmp33 + tmp36;
     }
     {
	  fftw_real tmp27;
	  fftw_real tmp28;
	  fftw_real tmp13;
	  fftw_real tmp24;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp27 = tmp25 + tmp26;
	  tmp28 = tmp18 + tmp23;
	  real_output[2 * real_ostride] = tmp27 - tmp28;
	  real_output[4 * real_ostride] = tmp27 + tmp28;
	  tmp13 = K866025403 * (tmp11 - tmp12);
	  tmp24 = tmp18 - tmp23;
	  imag_output[imag_ostride] = tmp13 - tmp24;
	  imag_output[5 * imag_ostride] = -(tmp13 + tmp24);
     }
}

fftw_codelet_desc fftw_real2hc_12_desc =
{
     "fftw_real2hc_12",
     (void (*)()) fftw_real2hc_12,
     12,
     FFTW_FORWARD,
     FFTW_REAL2HC,
     266,
     0,
     (const int *) 0,
};

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt