/* Submodel PNVO001 skeleton created by AME Submodel editing utility
   ven. 6. oct. 11:10:58 2017 */



#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "ameutils.h"
/* *******************************************************************************
TITLE : PNVO001
------------------------------------------------------------------------------
DATE OF CREATION / AUTHOR : 
   2002 : Created by FS from PNVO01
------------------------------------------------------------------------------
SOURCE :

   This material contains trade secrets or otherwise confidential
   information owned by Siemens Industry Software Inc. or its
   affiliates (collectively, "Siemens"), or its licensors. Access to
   and use of this information is strictly limited as set forth in the
   Customer's applicable agreements with Siemens.
   
   Unpublished work. Copyright 2023 Siemens

******************************************************************************* */

#define _SUBMODELNAME_ "PNVO001"

/* >>>>>>>>>>>>Insert Private Code Here. */

/* real stores */
#define PATM      0
#define AREAMAX   1
#define CQ        2

/* integer stores */
#define DISC_LIMIT   0
#define DISC_ORIF    1

/* <<<<<<<<<<<<End of Private Code. */


/* There are 4 real parameters:

   cq    flow coefficient (Cq)           [null]
   area0 orifice area at maximum opening [mm**2 -> m**2]
   Cv    maximum flow coefficient (Cv)   [null]
   Kv    maximum flow coefficient (Kv)   [null]
*/


/* There are 2 integer parameters:

   gi      gas type index          
   flowset flow coefficient setting
*/

void pnvo001in_(int *n, double rp[4], int ip[2], double c[3]
      , int ic[2])

{
   int loop, error;
/* >>>>>>>>>>>>Extra Initialization Function Declarations Here. */
/* <<<<<<<<<<<<End of Extra Initialization declarations. */
   int gi, flowset;
   double cq, area0, Cv, Kv;

   gi         = ip[0];
   flowset    = ip[1];

   cq         = rp[0];
   area0      = rp[1];
   Cv         = rp[2];
   Kv         = rp[3];
   loop = 0;
   error = 0;

/*
   If necessary, check values of the following:

   rp[0..3]
*/


/* >>>>>>>>>>>>Initialization Function Check Statements. */

   pn2_valid_gas_(&gi, &error);

   if (flowset == 1)
   {
      if (area0 < 0.0)
      {
         error = 2;
         amefprintf(stderr, "\nOrifice area at maximum opening should be positive.\n");
      }

      if (cq <= 0.0 )
      {
         error = 2;
         amefprintf(stderr, "\nFlow coefficient should be strictly positive.\n");
      }
   }
   else if (flowset == 2)
   {
      if (Cv < 0.0)
      {
         error = 2;
         amefprintf(stderr, "\nMaximum flow coefficient (Cv) should be positive (value is %g).\n", Cv);
      }
   }
   else
   {
      if (Kv < 0.0)
      {
         error = 2;
         amefprintf(stderr, "\nMaximum flow coefficient (Kv) should be positive (value is %g).\n", Kv);
      }
   }

/* <<<<<<<<<<<<End of Initialization Check Statements. */

/*   Integer parameter checking:   */

   if (gi < 1 || gi > 99)
   {
      amefprintf(stderr, "\ngas type index must be in range [1..99].\n");
      error = 2;
   }
   if (flowset < 1 || flowset > 3)
   {
      amefprintf(stderr, "\nflow coefficient setting must be in range [1..3].\n");
      error = 2;
   }

   SUBMODEL_HANDLE_AND_RESET_ERROR(_SUBMODELNAME_, n, error)

/* Common -> SI units conversions. */

   rp[1]    *= 1.00000000000000e-006;
   area0      = rp[1];


/* >>>>>>>>>>>>Initialization Function Executable Statements. */

   /* get atmospheric pressure */
   c[PATM] = pn2getatp_();

   if (flowset == 1)
   {
      c[CQ] = cq;
      c[AREAMAX] = area0;
   }
   else
   {
      /* calculation of equivalent maximal area with Cv or Kv.
         Default value of cq; the same value will be used in pn2rcqfix. */
      c[CQ] = 0.72;

      if (flowset == 2)   /* Cv */
         orif_areafromcv_(&Cv, &c[CQ], &c[AREAMAX]);
      else
         orif_areafromkv_(&Kv, &c[CQ], &c[AREAMAX]);
   }

/* <<<<<<<<<<<<End of Initialization Executable Statements. */
}

/*  There are 3 ports.

   Port 1 has 1 variable:

      1 res     input signal [null] basic variable input

   Port 2 has 4 variables:

      1 dh2        enthalpy flow rate at port 2 [J/s -> W]    basic variable output
      2 dm2        mass flow rate at port 2     [g/s -> kg/s] basic variable output
      3 temp2      temperature at port 2        [K]           basic variable input
      4 press2     pressure at port 2           [Pa]          basic variable input

   Port 3 has 4 variables:

      1 dh3        duplicate of dh2             (sign reversed)
      2 dm3        duplicate of dm2             (sign reversed)
      3 temp3      temperature at port 3        [K]  basic variable input
      4 press3     pressure at port 3           [Pa] basic variable input
*/

/*  There are 3 internal variables.

      1 xv         fractional opening          [null]            basic variable
      2 cm         mass flow parameter (cm)    [(kg*K/J)**(1/2)] basic variable
      3 gasvel     vena contracta gas velocity [m/s]             basic variable
*/

void pnvo001_(int *n, double *res, double *dh2, double *dm2
      , double *temp2, double *press2, double *temp3, double *press3
      , double *xv, double *cm, double *gasvel, double rp[4]
      , int ip[2], double c[3], int ic[2])

{
   int loop;
/* >>>>>>>>>>>>Extra Calculation Function Declarations Here. */
   double marea;   /* modulated area */
   double pressa2, pressa3;
   static double zero = 0.0, one = 1.0;
/* <<<<<<<<<<<<End of Extra Calculation declarations. */
   int gi, flowset;
   double cq, area0, Cv, Kv;

   gi         = ip[0];
   flowset    = ip[1];

   cq         = rp[0];
   area0      = rp[1];
   Cv         = rp[2];
   Kv         = rp[3];
   loop = 0;

/*
   Set all submodel outputs below:

   *dh2        = ??;
   *dm2        = ??;
   *xv         = ??;
   *cm         = ??;
   *gasvel     = ??;
*/



/* >>>>>>>>>>>>Calculation Function Executable Statements. */

   /* set absolute pressure */
   pressa2 = *press2 + c[PATM];
   pressa3 = *press3 + c[PATM];

   *xv = dlimit_(res, &zero, &one, &ic[DISC_LIMIT]);

   /* limitation of the modulated area */
   marea = *xv * c[AREAMAX];

   /*** calculation of the flows ***/
   pn2rcqfix_( dh2, dm2, temp2, &pressa2, temp3, &pressa3, &marea, &c[CQ], &gi,
               cm, gasvel, &ic[DISC_ORIF] );

/* <<<<<<<<<<<<End of Calculation Executable Statements. */

/* SI -> Common units conversions. */

   *dm2      /= 1.00000000000000e-003;
}